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\  ABSTRACT 

Vortex  aerodynamics  has  played  an  important  role  in  the  de¬ 
velopment  of  high  performance  aircraft  in  recent  years.  Although 
computer  codes  which  solve  the  three  dimensional  Euler  equations 
have  been  used  extensively  to  study  leading-edge  vortices,  they 
don’t  include  physical  viscosity  effects  associated  with  vortex 
flows.  The  Euler  solvers  do,  however,  contain  numerical  viscosi¬ 
ty.  As  a  result,  viscosity  effects  in  the  Euler  solutions  such 
as  vortex  core  size,  vortex  burst  location,  leading  edge  separa¬ 
tion,  and  vortex  rollup  often  do  not  agree  quantitatively  with 
results  of  physical (experiments.  The  present  work  defines  models 
for  these  phys i ca 1 ' v i scos i ty  effects  which  can  be  coupled  with  an 
Euler  solver  to  improve  modeling  of  vortex  physics. 

A  vortex  core  model  is  derived  from  the  steady,  incompress¬ 
ible  Naviei — Stokes  equations  written  in  cylindrical  coordinates. 

The  core  model  is  coupled  with  an  Euler  solver  and  tested  on  a 
variety  of  delta  wings  over  a  range  of  angles  of  attack.  The 
resulting  surface  pressure  distributions  and  vortex  burst  loca¬ 
tions  are  shown  to  be  much  closer  to  wind  tunnel  data  and  results 
from  Nav i er-Stokes  solutions  than  results  from  Euler  codes  alone,  U 

A  second  model  is  defined  for  viscosity  effects  in  the  vis-  -  ; 

cous  shear  layer  near  the  rounded  leading  edge  of  a  highly  swept  " 

wing  based  on  an  analogy  to  the  boundary  layer  on  a  flat  plate. 

The  model  is  incorporated  into  an  Euler  code  through  the  surface  .  •. 
boundary  condition  and  source  terms  in  the  boundary  cells.  The 
modified  code  is  tested  on  several  highly  swept  wings  with  round-' 
ed  leading  edges.  Results  are  also  shown  to  be  in  closer  agree¬ 
ment  with  wind  tunnel  data  for  the  same  wing  geometry  than  re- 
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ABSTRACT 


Vortex  aerodynamics  has  played  an  important  role  in  the  development  of 
high  performance  aircraft  in  recent  years.  Although  computer  codes  which 
solve  the  three  dimensional  Euler  equations  have  been  used  extensively  to 
study  vortex  flows,  they  don't  include  physical  viscosity  effects 
associated  with  vortex  flows.  The  Euler  solvers  do,  however,  contain 
numerical  viscosity.  As  a  result,  viscosity  effects  in  the  Euler  solutions 
such  as  vortex  core  size,  vortex  burst  location,  leading  edge  separation, 
and  vortex  rollup  often  do  not  agree  quantitatively  with  results  of 
physical  experiments.  The  present  work  defines  models  for  these  physical 
viscosity  effects  which  can  be  coupled  with  an  Euler  solver  to  improve 
modeling  of  vortex  physics. 

A  vortex  core  model  is  derived  from  the  steady,  incompressible 
Navier-Stokes  equations  written  in  cylindrical  coordinates.  The  core 
model  is  coupled  with  an  Euler  solver  and  tested  on  a  variety  of  delta 
wings  over  a  range  of  angles  of  attack.  The  resulting  surface  pressure 
distributions  and  vortex  burst  locations  are  shown  to  be  much  closer  to 
wind  tunnel  data  and  results  from  Navier-Stokes  solutions  than  results  from 
Euler  codes  alone . 

A  second  model  is  defined  for  viscosity  effects  in  the  viscous  shear 
layer  near  the  rounded  leading  edge  of  a  highly  swept  wing  based  on  an 
analogy  to  the  boundary  layer  on  a  flat  plate .  The  model  is  incorporated 
into  an  Euler  code  through  the  surface  boundary  condition  and  source  terms 
in  the  boundary  cells.  The  modified  code  is  tested  on  several  highly  swept 
wings  with  rounded  leading  edges.  Results  are  also  shown  to  be  in  closer 
agreement  with  wind  tunnel  data  for  the  same  wing  geometry  than  results 
from  an  unmodified  Euler  code. 
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CHAPTER  1 


INTRODUCTION 


1.1  MOTIVATION  OF  THE  PRESENT  WORK 


Vortex  aerodynamics  has  played  an  increasingly  important  role  in  the 
design  of  high  performance  aircraft  in  recent  years.  In  the  late  1940' s, 
wind  tunnel  studies  of  wings  with  triangular  or  delta  planforms  revealed 
that  such  wings  produced  more  lift  at  high  angles  of  attack  than  was 
predicted  by  conventional  finite  wing  theory.  The  first  of  these  studies, 
which  was  reported  by  Wilson  and  Lovell  [1],  found  that  the  magnitude  of 
this  extra  lift  was  increased  when  the  test  wing  was  given  a  sharp  leading 
edge.  Flow  visualization  in  wind  tunnels  and  water  channels  revealed  that 
as  a  delta  wing  approaches  an  angle  of  attack  where  flow  separation  would 
cause  an  unswept  wing  to  stall,  the  flow  over  each  highly  swept  leading 
edge  instead  rolls  up  into  a  tornado-like  flow  structure  called  a  leading 
edge  vortex.  The  use  of  a  sharp  leading  edge  causes  the  flow  to  separate 
and  form  a  vortex  at  a  lower  angle  of  attack.  Since  the  separated  flow 
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rolls  up  instead  of  being  disorganized  into  turbulence,  the  loss  of  lift 
normally  associated  with  a  stall  is  delayed  until  higher  angles  of  attack. 
Like  a  tornado,  leading  edge  vortices  have  very  low  pressures  in  their 
center  or  core  region.  Because  of  this,  they  lower  the  pressure  on  a 
portion  of  the  upper  surface  of  the  wing  and  hence  increase  the  amount  of 
lift  produced.  The  extra  vortex  lift  gives  delta-winged  aircraft  like  the 
F-102  and  F-106  extremely  good  maneuverability  at  low  speeds. 

The  benefits  of  leading  edge  vortices  have  prompted  designers  to 
include  vortex  lift-generating  surfaces  in  such  modern  high  performance 
aircraft  as  the  F-16,  F-18,  SR- 71,  Concorde  SST,  and  even  the  Space 
Shuttle.  Figure  1  shows  planform  views  of  these  aircraft  with  the 
leading-edge  vortex- generating  surfaces  identified  by  crosshatching.  Such 
aircraft  are  able  to  fly  at  angles  of  attack  greater  than  30  degrees  in 
some  cases  without  experiencing  a  stall.  Aircraft  with  better  high  speed 
performance  can  be  designed  by  relying  on  vortex  lift  for  their  low  speed, 
high  angle  of  attack  operations.  However,  care  must  be  taken  when 
designing  such  aircraft  to  avoid  problems  which  can  result  from  the 
presence  of  vortices  in  the  flow  field. 

As  more  experimental  data  on  leading  edge  vortices  were  accumulated, 
researchers  such  as  Hummel  [2]  and  Anders  [3]  found  that  when  a  delta  wing 
operates  at  very  high  angles  of  attack,  axial  pressure  gradients  develop  in 
the  cores  of  its  leading  edge  vortices.  As  angle  of  attack  is  increased, 
these  pressure  gradients  eventually  cause  local  stagnation  or  reversal  of 
the  flow  in  the  core.  When  this  happens,  it  causes  the  core  to  grow 
rapidly  and  break  down  into  random  turbulence  in  much  the  same  way  that 
stall  occurs  on  unswept  wings.  This  phenomenon  is  called  vortex  bursting. 
Several  researchers  including  Wentz  [4],  Payne  and  Nelson  [5],  Hall  [6], 
and  Lawford  and  Beauchamp  [7]  have  made  extensive  wind  tunnel  and  water 
channel  tests  of  a  variety  of  delta  wings.  They  have  created  a  large  body 
of  bursting  location  data  as  a  function  of  sweep  angle  and  angle  of  attack. 
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Figure  1.  Planfora  views  of  some  high  performance  aircraft.  Surfaces 
generate  leading  edge  vortices  are  crosshatched. 
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When  bursting  occurs,  lift  is  greatly  reduced.  If  the  vortex  on  one 
wing  of  an  aircraft  bursts  before  that  on  the  other  wing,  the  aircraft 
experiences  severe  control  problems.  Because  bursting  progresses  from  the 
wing  trailing  edge  toward  the  apex,  it  can  also  cause  an  uncontrollable 
pitchup  on  some  delta  wings.  In  addition,  the  turbulence  present  in  a 
burst  vortex  can  cause  structural  fatigue  problems  in  aircraft  components 
which  are  not  designed  for  that  frequency  of  cyclic  loading.  Skow, 

Titriga,  and  Moore  [8],  among  others,  have  described  aircraft 
configurations  which  even  experience  control  difficulties  due  to  the  motion 
and  interaction  of  the  vortices  they  produce  before  any  bursting  occurs. 

Present  technologies  for  understanding  vortices  and  utilizing  them  in 
aircraft  design  have  some  serious  limitations.  Wind  tunnels  and  water 
channels  are  currently  used  for  the  majority  of  vortex-related  design 
development  work,  but  such  testing  is  expensive  and  requires  long  lead 
times,  making  it  poorly  suited  to  preliminary,  iterative  design  work.  Wind 
tunnels  also  provide  only  limited  access  to  local  flow  details. 

Computational  fluid  dynamics  (CFD)  technologies  have  shown  consider¬ 
able  ability  to  model  vortex  effects.  Changes  to  the  flow  conditions  and 
vehicle  configurations  being  modeled  are  accomplished  in  CFD  codes  with 
less  effort  than  is  required  for  changes  in  wind  tunnel  models  and  test 
conditions.  This  makes  CFD  very  attractive  for  preliminary  design  and 
configuration  definition  work. 

Several  mathematical  models  and  numerical  methods  have  been  developed 
to  simulate  leading  edge  vortices.  As  available  computing  power  has 
increased,  solution  technology  has  progressed  from  the  leading  edge  suction 
analogy  of  Polhamus  [9]  to  panel  methods,  [10,11]  then  to  Euler  equation 
formulations.  [12,13,14,15]  With  each  new  method,  the  requirement  for 
computer  speed  and  memory  has  increased.  In  recent  years,  attempts  have 
been  made  to  solve  the  Navier-Stokes  equations  for  flowfields  around  delta 


wings.  [16,17,18]  This  technology,  however,  is  not  yet  mature.  Numerical 
problems  associated  with  grid  density  and  truncation  error  where  high 
gradients  are  present  have  yet  to  be  resolved.  In  addition,  current 
Navier- Stokes  codes  employ  turbulence  models  which  are  inadequate  for 
vortex  flows.  Computing  time  and  memory  requirements  of  Navier-Stokes 
codes  also  make  them  too  expensive  for  routine  uses. 

On  the  other  hand,  computer  codes  which  solve  the  Euler  equations 
offer  a  low-cost  alternative  to  the  Navier-Stokes  solvers  for  problems 
where  viscosity  effects  are  either  negligible  or  can  be  represented  by 
simple  mathematical  models.  Many  researchers,  including  Rizzi  [19]  and 
Newsome  [20] ,  have  made  extensive  use  of  Euler  solvers  to  simulate  leading 
edge  vortex  flows.  Their  results  are  often  in  good  qualitative  agreement 
with  experimental  and  Navier-Stokes  results,  but  they  frequently  fail  to 
accurately  predict  such  viscosity  effects  as  vortex  core  size  and  bursting 
location.  Analytic  models  for  physical  viscosity  effects  can  be  coupled 
with  the  Euler  codes.  The  formulation  of  such  a  method  would  contribute  to 
a  better  understanding  of  the  processes  which  occur  in  a  leading  edge 
vortex.  The  combination  of  Euler  code  and  model  would  be  a  useful  tool  for 
vortex  research  and  aircraft  design.  This  is  the  motivation  for  the 
present  research. 


1.2  DESCRIPTION  OF  THE  PROBLEM 


Computer  codes  which  solve  the  Euler  equations  do  not  attempt  to 
include  physical  viscosity  effects.  Ideally,  an  Euler  solution  should  be 
totally  inviscid.  In  practice,  however,  numerical  effects  which  are 
similar  to  physical  viscosity  are  present  in  Euler  solutions.  This 
numerical  viscosity  will  be  described  in  greater  detail  in  Chapter  2. 

Because  of  numerical  viscosity,  Euler  solutions  contain  such  viscous 
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phenomena  as  flow  separation,  leading  edge  vortices,  and  vortex  breakdown. 
These  effects  are  not,  however,  controlled  by  the  viscous  physics  of  the 
flow  being  modeled.  For  flow  problems  where  viscosity  effects  are 
relatively  unimportant,  the  difference  between  numerical  and  physical 
viscosity  has  a  relatively  small  effect  on  the  flow  solutions.  For  these 
cases,  Euler  solvers  often  give  acceptable  results. 

In  flow  problems  involving  leading  edge  vortices,  however,  viscosity 
effects  can  have  a  profound  influence  on  the  very  structure  of  the  flow 
solution.  In  particular,  the  leading  edge  separation  which  leads  to  vortex 
rollup  and  the  onset  of  vortex  bursting  are  both  strongly  affected  by 
viscosity.  In  general,  numerical  viscosity  effects  are  in  the  same 
direction  as  physical  viscosity,  so  there  is  often  qualitative  agreement 
between  Euler  solutions  and  physical  experiments.  The  magnitudes  of 
physical  and  numerical  viscosity  frequently  differ,  however.  While 
physical  viscosity  is  solely  a  function  of  flow  quantities,  numerical 
viscosity  is  also  influenced  by  the  computational  algorithm  used  and  the 
size,  shape  and  orientation  of  the  grid.  As  a  result,  quantitative 
agreement  between  computations  and  physics  is  often  difficult  to  achieve. 

A  need  therefore  exists  for  a  method  to  replace  or  correct  the  effects 
of  numerical  viscosity  in  Euler  solutions  so  that  they  better  represent  the 
viscous  physics.  In  the  present  effort,  viscosity  effects  are  measured  by 
flow  quantities  in  the  Euler  solution  such  as  pressure,  velocity,  vortex 
strength  and  core  diameter.  Corrections  to  the  flow  quantities  are  then 
calculated  from  algebraic  models  for  the  viscous  processes  near  the  wing 
surface  and  in  the  vortex  core.  These  corrections  are  incorporated  into 
the  Euler  code  as  source  terms  in  the  flow  equations.  For  this  method  to 
succeed,  the  algebraic  models  must  properly  represent  the  relevant  viscous 
physics . 
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1.3  EXISTING  VORTEX  MODELS 


A  great  deal  of  theoretical  modeling  of  vortex  decay  and  vortex 
bursting  has  been  done  by  a  variety  of  researchers  in  the  past  50  years. 
Research  on  this  topic  has  typically  been  divided  into  two  areas  of 
interest.  The  oldest  studies  of  vortex  decay  dealt  with  aircraft  trailing 
vortices.  When  vortex  bursting  on  delta  wings  was  discovered  in  the  late 
40' s,  a  great  deal  of  attention  was  given  to  trying  to  understand  this 
phenomenon.  In  the  late  60' s  and  early  70' s,  when  jumbo  jets  were 
introduced  into  commercial  airline  service,  their  trailing  vortices  created 
serious  hazards  to  other  aircraft.  Attention  once  more  turned  to  the 
modeling  of  trailing  vortex  decay.  Models  for  vortex  bursting  were  studied 
in  search  of  a  means  to  alleviate  the  hazard  of  trailing  vortices  by 
causing  them  to  burst.  Then,  in  the  late  70's  and  80's  the  design  of  high 
performance  aircraft  such  as  the  F-16,  F-18,  and  Space  Shuttle  focused 
interest  once  more  on  leading  edge  vortex  bursting. 

All  of  the  models  and  theories  developed  in  this  area  have 
implications  for  the  present  research.  Each  model  has  certain  advantages 
and  disadvantages  for  describing  vortex  bursting.  A  survey  of  all  the 
models  reveals  at  the  same  time  a  great  variety  in  their  formulation  and 
yet  a  general  similarity  in  their  fundamental  assumptions ,  especially  the 
shapes  of  their  assumed  velocity  distributions. 

One  of  the  earliest  and  simplest  3-dimensional  vortex  models  was  the 
potential  vortex.  Circumferential  velocity  in  this  vortex  varies  with  the 
inverse  of  the  distance  from  the  center.  This  causes  circulation  around 
all  closed  paths  containing  the  vortex  center  to  be  the  same.  The 
potential  vortex  model  has  no  viscous  core,  so  circumferential  velocity 
becomes  infinite  at  the  center.  Far  from  the  center,  however,  the  velocity 
field  matches  well  with  those  of  vortices  in  nature.  Most  subsequent 
vortex  models  have  been  designed  to  match  the  potential  vortex  velocity 
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field  outside  their  core  regions.  Axial  velocities  for  the  potential 
vortex  are  defined  as  being  uniform. 

A  simple  modification  of  the  potential  vortex  avoids  the  singularity 
at  the  center.  Known  as  the  Rankine  vortex,  this  model  has  a  central  core 
which  rotates  as  a  solid  body.  At  the  core  edge,  the  circumferential 
velocity  profile  abruptly  changes  to  that  of  the  potential  vortex. 
Circumferential  velocities  match  at  the  core  edge,  but  slopes  do  not.  The 
great  advantage  of  the  Rankine  vortex  for  modeling  is  the  simplicity  of 
this  circumferential  profile.  For  this  reason,  it  has  been  a  popular 
velocity  profile  for  the  development  of  bursting  theories  and  simulations. 
Like  the  potential  vortex,  axial  velocities  for  the  Rankine  vortex  are 
uniform. 

A  more  physically  reasonable  vortex  model  was  derived  by  Lamb  [23]  as 
a  similarity  solution  for  the  viscous  decay  of  an  infinite  line  vortex 
which  is  initially  concentrated  into  a  potential  vortex.  The  circumferen¬ 
tial  velocities  in  this  model  are  expressed  as  a  function  of  time.  For 
this  reason,  the  Lamb  vortex  was  used  by  Iversen  [24]  to  model  trailing 
vortex  decay.  In  Iversen' s  model,  a  variable  turbulent  eddy  viscosity  is 
used  instead  of  Lamb's  constant  viscosity.  This  allows  Iversen  to  match 
water  channel  trailing  vortex  decay  data  more  accurately.  Both  models  do 
not  consider  axial  variations  in  vortex  properties.  This  prevents  them 
from  modeling  vortex  bursting  in  their  present  form. 

The  analytic  modeling  of  vortex  bursting  really  started  with  Hall's 
[25]  work  published  in  1961.  He  assumed  that  once  viscosity  effects  had 
determined  the  circumferential  velocity  profile  of  the  vortex,  the  bursting 
process  could  by  modeled  inviscidly.  Hall  also  assumed  a  conical, 
incompressible  vortex  with  steady,  laminar  flow.  The  initial  velocity 
profiles  assumed  by  Hall  allowed  variation  in  the  axial  velocities,  but 
both  his  axial  and  circumferential  profiles  had  singularities  at  the  vortex 
center.  This  forced  him  to  define  a  diffusive  sub-core  to  avoid  the 
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singularity  and  then  patch  the  two  solutions.  As  a  result,  Hall's  model 
was  relatively  complex. 

Other  researchers  such  as  Brown,  [26]  Mager,  [27]  Krause  [28]  and 
Powell  and  Murman  [29]  have  derived  similar  sets  of  differential  equations 
to  describe  bursting.  Each  used  slightly  different  assumptions  and 
methods,  but  all  used  some  arbitrarily  defined  initial  velocity 
distribution.  The  various  velocity  profiles  employed  are  described  and 
compared  in  Appendix  A.  Each  researcher  solved  the  set  of  equations 
numerically  and/or  analytically.  Mager  obtained  a  closed  form  solution  for 
his  ordinary  differential  equations  for  the  case  of  an  isolated  vortex  in  a 
uniform  flow.  Because  the  flow  field  around  a  delta  wing  is  far  from 
uniform,  Mager 's  closed  form  solution  is  not  used  in  the  present  work. 

Both  Krause  [30]  and  Grabowski  and  Berger  [31]  used  the  algebraic  velocity 
profiles  defined  by  Mager  to  make  numerical  simulations  of  vortex  bursting. 
The  results  of  these  simulations  agree  well  with  bursting  experiments 
performed  in  vortex  tubes  by  Harvey  [32],  Sarpkaya  [33],  and  Faler  and 
Leibovich  [34], 

A  great  deal  of  theoretical  and  experimental  work  has  been  done  aimed 
at  explaining  vortex  bursting  in  terms  of  stability.  The  first  of  these 
was  Jones  [34],  who  analyzed  the  response  of  an  arbitrary  cylindrical  flow 
to  infinitesimal  axial  disturbances.  Jones  then  applied  his  results  to 
Hall's  solution.  Subsequent  analyses  by  Benjamin  [35],  Lessen  [36], 
Leibovich  [37]  and  others  have  described  vortex  breakdown  as  a  transition 
or  jump  between  two  flow  states,  one  of  which  has  reverse  axial  flow,  such 
that  both  flow  states  satisfy  the  governing  differential  equations.  They 
defined  stability  bounds  which  suggest  that  some  vortices  with  "super¬ 
critical"  initial  axial  and  circumferential  velocity  distributions  are 
prone  to  bursting  while  vortices  with  "subcritical"  initial  profiles  are 
not.  The  numerical  simulations  by  Grabowski  and  Berger  tended  to  disagree 
with  these  stability  bounds,  however.  Mager  discovered  similar  conjugate 
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flow  states  in  his  closed  form  solutions.  He  found  that  the  change  from 
one  flow  state  to  another  could  either  be  a  smooth  transition  or  an  abrupt 
jump,  depending  on  the  nature  of  the  initial  velocity  profile. 

All  of  these  methods  have  used  arbitrary  initial  velocity  profiles  as 
part  of  their  analyses.  These  velocity  profiles  have  generally  been 
axisymmetric .  The  flow  external  to  the  vortex  core  has  been  modeled  as 
either  a  uniform  velocity  field  with  a  line  vortex  added  to  it,  or  as  a 
conical  flow.  For  these  reasons,  strict  applicability  of  any  of  the 
solutions  or  criteria  to  a  general  leading  edge  vortex  problem  is  limited. 
Some  of  the  basic  methods  used  to  obtain  these  solutions  are,  however,  used 
in  the  present  work. 


1.4  OUTLINE  OF  THE  PRESENT  WORK 


A  model  for  viscosity  effects  in  vortex  cores  is  derived  from  the 
steady,  incompressible  Navier-Stokes  equations  written  in  cylindrical 
coordinates.  The  model  is  formulated  to  allow  it  to  be  incorporated  into  a 
computer  code  which  solves  the  Euler  equations  in  three  dimensions.  The 
resulting  Euler  code/core  model  combination  is  tested  on  a  variety  of 
problems.  The  first  tests  are  made  on  a  simple  rectangular  box  domain. 
These  tests  determine  the  ability  of  the  core  model  to  correct  for  the 
effects  of  numerical  viscosity  in  the  Euler  solution.  Further  tests  are 
made  on  the  rectangular  grid  to  determine  the  model's  ability  to  produce 
vortex  bursting  when  a  rapidly  growing  physical  core  is  modeled  and  to 
prevent  bursting  if  a  physical  core  with  zero  growth  is  assumed.  The  Euler 
code  with  core  model  is  then  tested  on  a  variety  of  delta  wings  over  a 
range  of  angles  of  attack.  The  resulting  surface  pressure  distributions 
and  vortex  burst  location  predictions  are  compared  with  wind  tunnel  data 
and  results  from  Navier-Stokes  simulations  and  Euler  codes  alone. 
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A  model  for  total  pressure  losses  due  to  viscosity  in  the  shear  layer 
near  the  rounded  leading  edge  of  a  highly  swept  wing  is  defined  based  on  an 
analogy  to  the  boundary  layer  on  a  flat  plate.  The  model  is  incorporated 
into  an  Euler  code  through  the  surface  boundary  condition  and  momentum 
source  terms  in  the  boundary  cells.  The  modified  code  is  tested  over  a 
range  of  angles  of  attack  for  several  highly  swept  wings  with  rounded 
leading  edges.  These  flow  solutions  are  compared  with  wind  tunnel  data  and 
results  from  an  unmodified  Euler  code. 
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CHAPTER  2 

COMPUTATIONAL  METHODS 


2.1  EULER  EQUATIONS 

The  Navier- Stokes  equations  are  generally  used  to  describe  the 
unsteady  motion  of  a  compressible  viscous  fluid.  For  many  flow  problems, 
however,  viscosity  effects  are  confined  to  relatively  small  portions  of  the 
flowfield.  For  instance,  for  low  speed  flow  around  an  infinite  wing, 
viscosity  effects  may  be  confined  to  a  thin  boundary  layer  close  to  the 
wing  surface.  If  this  region  is  modeled  appropriately,  the  rest  of  the 
flowfield  may  be  solved  with  a  simplified  version  of  the  Navier-Stokes 
equations  in  which  the  viscosity  terms  are  omitted.  These  simplified 
equations  are  known  as  the  Euler  equations.  The  use  of  the  Euler  equations 
bypasses  the  computational  difficulties  of  dealing  with  the  very  small 
length  scales  associated  with  viscosity  effects. 

A  number  of  computer  codes  are  available  to  solve  the  Euler  equations 
for  flowfields  around  three-dimensional  aerodynamic  shapes.  Two  of  these 
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codes,  FL057  [13]  and  ARC3D  [15],  are  modified  and  used  in  the  present 
research.  The  computational  methods  employed  in  each  of  these  codes  will 
be  described  because  the  present  research  is  directed  at  overcoming  some  of 
the  limitations  inherent  in  those  methods. 

The  time  dependent,  compressible  Euler  equations  can  be  written  in 
matrix  form: 


dui  dj  dg  dh 
~dt  +  Tx  +  dY  +  ~dZ 
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and  U,  V,  and  W  are  component  velocities  in  the  three  orthogonal  coordinate 
directions,  X,  Y,  and  Z  of  a  Cartesian  coordinate  system.  Also,  p  is  the 
local  air  density,  p  is  the  local  air  pressure,  and  E  and  H  are  given  by: 
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2.2  FLO 5 7  -  EXPLICIT  FINITE  VOLUME  FORMULATION 


In  FL057  [13],  the  time  dependent  Euler  equations  are  solved  to  obtain 
steady-state  solutions  as  time  asymptotes.  The  discretization  scheme  used 
by  FL057  is  known  as  the  finite  volume  method.  In  this  method,  the  flow 
field  is  divided  into  small  cells  or  control  volumes.  The  Euler  equations 
are  satisfied  in  the  average  sense  for  each  control  volume.  This  is 
essentially  solving  the  integral  form  of  (1)  which  is 


/ 


/ 


dui 

Jt 


dj  dg  dh 
dX+dY+dZ 


dX  dY  dZ  =  0  . 


(4) 


and  which,  for  small  cells,  can  be  rewritten: 


—  (wh)  +  Qui=  0 
dt 


(5) 


where  h  is  the  cell  volume  and  Q  is  the  operator  for  the  steady- state 
residual  calculation  applied  to  the  flow  quantities.  In  other  words,  the 
rate  of  change  of  the  flow  quantities  in  each  cell  is  equated  to  the  net 
fluxes  of  those  quantities  across  the  cell  walls.  For  the  finite  volume 
method,  all  of  the  flow  quantities  are  defined  at  the  center  of  the  cells, 
rather  than  at  the  grid  nodes  at  the  cell  corners.  When  evaluating  fluxes, 
however,  flow  quantities  on  each  face  of  the  cell  must  be  known.  These 
quantities  are  approximated  by  making  their  values  at  the  cell  face  equal 
to  the  average  of  quantities  at  the  centers  of  the  cells  which  share  that 
face.  Derivatives  of  flow  quantities  are  approximated  as  differences 
between  the  quantities  in  adjacent  cells. 

In  order  to  stabilize  the  solution,  numerical  viscosity  is  added  to 
the  formulation.  In  FL057,  the  numerical  viscosity  is  scaled  by  the 
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magnitudes  of  blended  second  and  fourth  derivatives  of  the  flow  quantities. 
With  this  damping  added,  the  final  form  of  (5)  which  is  actually  solved 
becomes : 


—  (wh.)  +  Qw-Dui*‘0 


(6) 


where  Dw  represents  the  sum  of  the  numerical  viscosity  added  for  stability 
plus  a  viscosity- like  numerical  effect  which  results  from  the 
discretization  of  the  flowfield  and  approximation  of  derivatives  with 
finite  differences.  For  simplicity  in  future  discussions,  the  term 
"numerical  viscosity"  will  be  used  in  this  work  to  refer  to  this  sum  of 
both  effects. 

The  flow  field  is  initialized  with  a  uniform  velocity  distribution, 
the  initial  condition,  and  then  a  fourth  order  Runge-Kutta  time-stepping 
scheme  is  used  to  advance  to  a  steady  state  solution.  The  Runge-Kutta 
scheme  uses  a  multi-level  explicit  integration  scheme  for  computational 
efficiency.  If  equation  (5)  is  rewritten 


dw  _  „ 
—  +  Pw 
dt 
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then,  at  time  level  n,  the  four  stage  scheme  can  be  written 
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A  single  iteration  involves  sweeping  through  the  entire  flowfield  and 
updating  each  flow  quantity.  The  Runge-Kutta  scheme  does  not  require  any 
matrix  inversion  for  integration,  but  it  allows  larger  time  steps  than 
other  explicit  Euler  codes.  Like  other  explicit  schemes,  the  method  also 
vectorizes  easily  on  supercomputers. 


2.3  ARC3D  -  IMPLICIT  FINITE  DIFFERENCE  FORMULATION 


Although  it  solves  the  same  Euler  equations  as  FL057,  the  details  of 
the  formulation  of  ARC3D  are  different.  The  major  features  which 
distinguish  ARC3D  from  FL057  are  a  finite  difference  formulation  instead  of 
finite  volume  and  an  implicit  time -stepping  scheme.  ARC3D  also  has  the 
capability  of  solving  the  thin-layer  approximation  to  the  Navier-Stokes 
equations  by  including  appropriate  viscous  terms. 

The  finite  difference  formulation  of  ARC3D  differs  from  finite  volume 
in  that  flow  quantities  are  defined  at  the  grid  node  points  instead  of  at 
the  cell  centers.  Like  FL057,  ARC3D  uses  a  centered  differencing  scheme  to 
obtain  spacial  derivatives.  Since  flow  quantities  are  not  weighted  by  cell 
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areas  as  in  FL057,  the  flow  equations  are  transformed  from  Cartesian 
coordinates  into  general  curvilinear  coordinates.  The  transformation 
metrics  are  generated  numerically  using  finite  differences. 

The  implicit  time -stepping  scheme  of  ARC3D  is  based  on  a  method 
developed  by  Beam  and  Warming  [22].  In  this  method,  the  equation  to  be 
solved  is 


d-*+rw~' 

dt 


0 


which  may  be  written 

uin"-uin  +  AtPti)n"  “  0. 

If  a  correction  is  defined  such  that 

6w  -  wn"-wn, 

then  the  equation  to  be  solved  can  then  be  written 

6ui+  AtP6w  =  (/  +  AtP)6u)  =  -  At  P  wn 

where  P  is  a  linearized  operator  of  P.  This  equation  can  be  written  as 
implicit  and  explicit  portions  to  yield 

N6w  =  -  At  Pw  (6) 

where  N  is  a  correction  operator  defined  in  three  dimensions  such  that 
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N 
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The  left  half  of  equation  (6)  is  implicit  while  the  right-hand  side  is  the 
explicit  evaluation  of  the  residual  at  time  step  n.  The  operator  N  is 
approximately  factored  into  three  one -dimensional  operators,  each  of  which 
can  be  inverted  using  a  block  tridiagonal  matrix  solver.  This  requires 
much  less  computational  effort  than  inverting  N,  but  computational  effort 
per  iteration  is  still  much  greater  than  in  explicit  schemes  like  FL057. 
This  extra  work  per  iteration  is  generally  offset  by  much  larger  time  steps 
permitted  by  the  implicit  scheme.  For  problems  considered  in  the  present 
work,  however,  ARC3D  typically  took  more  computer  time  than  FL057  to  reach 
the  same  level  of  convergence  in  the  flow  solution. 

As  in  FL057,  numerical  viscosity  is  added  based  on  second  and  fourth 
derivatives  of  flow  quantities.  The  second  order  dissipation  is  added 
implicitly  to  the  left  side  of  equation  (6)  to  improve  the  practical 
stability  bound.  The  fourth  order  term  is  added  explicitly  to  the 
formulation. 

The  thin- layer  approximations  to  the  Navier-Stokes  equations  involve 
neglecting  viscosity  terms  based  on  derivatives  in  the  streamwise 
direction.  As  a  result,  they  are  not  really  valid  for  the  case  of  a  vortex 
core.  Likewise,  separated  flow  invalidates  the  thin- layer  assumptions. 

For  these  reasons,  the  Navier-Stokes  capabilities  of  ARC3D  are  not  used  to 
any  great  extent  in  the  present  effort. 


2.4  BOUNDARY  CONDITIONS 


Appropriate  boundary  conditions  must  be  specified  to  make  a 
formulation  complete.  Three  different  types  of  boundary  conditions  are 
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implemented  in  the  present  work,  farfield,  solid  wall  and  symmetry.  At  far 
field  boundaries,  just  setting  the  flow  quantities  to  free  stream  values  is 
not  normally  acceptable.  If  such  an  approach  is  taken,  the  fixed 
quantities  at  the  boundaries  cause  wave  reflections  at  the  boundary  which 
propagate  back  into  the  flowfield.  This  can  lead  to  errors  and  instability 
in  the  solution. 

Instead  of  fixed  free  stream  conditions  at  farfield  boundaries,  a 
characteristic  boundary  condition  is  used  which  is  based  on  the  theory  of 
Kreiss  [21].  If  the  flow  is  subsonic,  this  theory  indicates  that  for 
boundaries  which  have  incoming  flow  through  them,  all  but  one  of  the  five 
needed  flow  quantities  must  be  specified  by  the  external  flow.  The  one 
remaining  quantity,  usually  pressure,  is  extrapolated  from  the  interior  so 
that  sound  waves  can  be  allowed  to  escape.  For  farfield  boundaries  with 
outflow,  on  the  other  hand,  only  one  quantity  is  specified  by  the  exterior 
flow  and  the  rest  are  extrapolated  from  the  interior. 

At  solid  surfaces,  flow  tangency  is  enforced  for  inviscid  flow 
problems.  This  is  done  by  setting  the  normal  component  of  the  flow 
velocity  at  the  surface  to  zero.  The  pressure  at  the  surface  must  also  be 
known  because  it  provides  a  momentum  source  to  the  boundary  cells  or  nodes . 
The  simplest  approximation  to  the  surface  pressure  is  to  equate  it  to  the 
pressure  at  the  boundary  cell  center  or  the  adjacent  node.  A  higher  order 
estimate  of  surface  pressure  can  be  achieved  by  using  the  momentum  equation 
normal  to  the  surface  to  extrapolate  pressure  from  the  cell  center  or  node 
down  to  the  solid  wall. 

For  planes  of  symmetry,  a  mirror  image  condition  for  all  the  flow 
quantities  is  enforced.  As  with  inviscid  solid  surfaces,  pressure  at  the 
plane  of  symmetry  may  be  extrapolated  from  the  interior  using  the  momentum 
equation  normal  to  the  symmetry  plane.  This  pressure  may  also  be  simply 
equated  to  the  pressure  at  the  adjacent  boundary  cell  center  or  grid  node. 
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2 . 5  GRIDS 


Three  types  of  grid  topology  are  used  in  the  present  research.  The 
first  of  these  is  a  simple  Cartesian  grid  in  a  rectangular  box  domain.  A 
perspective  view  of  one  such  grid  is  shown  in  Figure  2.  The  topology  of 
such  a  grid  is  labeled  H-H  because  boundaries  at  opposite  ends  of  the  box 
are  parallel  to  each  other  like  the  two  vertical  legs  of  the  letter  H. 

In  discussing  grid  topology,  it  is  helpful  to  refer  to  the  i,  j,  k 
indices  of  a  grid  point  instead  of  its  physical  coordinates,  X,  Y  and  Z. 

For  Cartesian  grids  with  unit  spacing  between  nodes,  there  is  no 
distinction  between  the  two  methods.  However,  as  a  grid  is  distorted  to 
conform  to  more  complex  geometries,  the  X,  Y,  and  Z  coordinates  of  a  given 
node  may  change,  but  the  computational  coordinates  given  by  its  i,  j  and  k 
indices  do  not.  The  convention  used  in  this  work  is  that  X  and  i 
directions  run  generally  in  the  streamwise  direction,  Y  and  j  run 
vertically,  and  Z  and  k  run  laterally  in  a  direction  that  makes  the  XYZ 
coordinate  system  right  handed.  In  FL057 ,  cells  are  indexed  in  the  same 
way  as  grid  nodes,  except  that  ghost  cells  are  added  along  the  boundaries 
of  the  computational  domain.  As  a  result,  the  dimensions  of  the  arrays 
storing  flow  quantities  are  larger  by  one  in  each  direction  than  those  for 
the  grid. 

Figure  3  shows  a  typical  boundary  cell  on  the  lower  surface  of  the 
rectangular  box.  Only  i  and  j  indices  are  shown  for  simplicity.  Note  that 
the  j  index  of  the  cell  is  2  because  the  ghost  cell  j  index  is  1.  The  grid 
nodes  at  the  surface  have  a  j  index  of  1  however.  Similar  relationships  at 
the  upstream  and  downstream  boundaries  cause  a  cell  with  a  given  i  index  to 
be  between  grid  nodes  with  indices  of  i-1  and  i.  The  H-H  grid  used  in  this 
work  has  an  inviscid  solid  surface  boundary  condition  on  the  j-1  face  of 
the  domain  and  a  symmetry  boundary  condition  on  the  k-1  face.  Farfield 
boundary  conditions  are  used  on  the  other  four  faces. 
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Figure  2 


Figure  3. 


Perspective  view  of  rectangular  box  grid. 
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A  second  grid  is  used  for  studies  of  leading  edge  vortices  on  delta 
wings.  In  this  grid,  the  i-1  and  i-imax  faces  remain  at  the  upstream  and 
downstream  ends  of  the  flowfield,  parallel  to  each  other  as  in  the  H-H 
grid.  However,  the  lower  surface  of  the  grid  is  wrapped  around  the  upper 
and  lower  surfaces  of  the  delta  wing,  causing  the  two  lateral  side  faces  of 
the  grid  at  k-1  and  k-kmax  to  end  up  both  being  on  the  plane  of  symmetry, 
one  on  the  top  side  of  the  wing  and  one  on  the  bottom.  In  the  region  of 
the  grid  ahead  of  the  apex  of  the  delta  wing,  the  j-1  surface  of  the  grid 
is  condensed  into  a  single  line.  A  perspective  view  of  this  grid  and  a 
planform  view  of  a  delta  wing  with  the  pattern  of  the  grid  on  its  surface 
are  shown  in  Figure  4. 

When  the  grid  on  the  other  side  of  the  plane  of  symmetry  is  included, 
the  j  -  constant  grid  lines  on  i  -  constant  planes  are  circles  or  ellipses. 
For  this  reason,  the  topology  of  this  grid  is  labeled  H-0.  The  H  implies 
that  i  -  constant  planes  are  still  parallel  as  in  the  Cartesian  grid  while 
the  0  suggests  the  circular  shape  of  the  j  -  constant  grid  lines. 

Boundary  conditions  for  this  grid  differ  from  those  of  the  H-H 
topology.  In  the  H-0  system,  both  k-1  and  k-kmax  faces  require  symmetry 
conditions.  The  j-1  surface  has  two  types  of  boundary  conditions.  For  i 
indices  running  from  the  apex  of  the  delta  wing  to  its  trailing  edge,  a 
solid  surface  boundary  condition  is  used.  Upstream  and  downstream  of  the 
wing,  a  periodic  boundary  condition  is  used  which  recognizes  the  fact  that 
a  cell  touching  the  j-1  surface  of  the  top  half  of  the  grid  is  physically 
adjacent  to  another  boundary  cell  with  a  higher  k  index  on  the  lower  half 
of  the  grid.  Figure  5  shows  a  view  of  an  i  -  constant  plane  of  the  grid 
downstream  of  the  wing  to  make  this  relationship  between  the  boundary  cells 
more  clear.  The  periodic  condition  can  be  handled  easily  by  setting  flow 
quantities  in  the  upper  cell  of  the  periodic  boundary  equal  to  the  flow 
quantities  in  the  corresponding  lower  cell  and  vice  versa  for  FL057,  and  by 
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Perspective  view  of  H-0  grid  and  planform  view  of  delta  wing. 
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Figure  5.  View  from  downstream  of  i  -  45  plane  of  H-0  grid. 
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making  an  appropriate  arrangement  between  the  upper  and  lower  grid  nodes 
for  ARC3D.  Farfield  boundary  conditions  are  used  for  the  rest  of  the 
boundaries . 

The  third  grid  topology  used  is  labeled  C-H.  As  the  name  implies,  the 
streamwise  direction  of  the  grid  is  folded  into  a  C  shape  while  the 
spanwise  direction  remains  similar  to  the  rectangular  box.  Figure  6  shows 
a  view  of  the  k-1  plane  along  with  a  wing  planform  view  and  grid  pattern 
for  a  C-H  grid  used  in  the  present  work.  The  C  shape  of  the  k-1  plane  is 
apparent.  It  is  also  apparent  from  Figure  6  that,  as  in  the  H-0  mesh,  two 
different  boundary  conditions  are  needed  for  the  j-1  surface. 

The  i  index  in  this  topology  starts  at  the  lower  half  of  the  down¬ 
stream  end  of  the  grid,  increases  moving  forward  and  wrapping  around  the 
wing  leading  edge,  and  ends  at  the  upper  half  of  the  downstream  end.  For 
those  cells  which  lie  on  the  wing  surface,  a  solid  surface  boundary 
condition  is  used.  A  periodic  boundary  condition  is  used  for  the  remaining 
cells  on  this  boundary.  This  includes  all  cells  on  this  surface  for  k 
greater  than  the  k  value  of  the  wingtip,  and  for  k  less  than  k  of  the 
wingtip,  all  cells  downstream  of  the  wing  trailing  edge.  Another  boundary 
condition  in  the  C-H  topology  which  differs  from  those  in  either  H-H  or  H-0 
meshes  is  the  i-1  boundary.  Since  the  i-1  surface  is  downstream  of  the 
wing  in  a  C-H  mesh,  it  is  given  an  outflow  farfield  boundary  condition. 


2.6  SPECIAL  CONCERNS  FOR  VORTEX  FLOWS 


When  a  vortex  is  present  in  an  Euler  solution,  numerical  viscosity  and 
grid  density  work  together  to  determine  vortex  core  size.  Rizzi(19)  made 
Euler  simulations  of  leading  edge  vortex  flows  on  delta  wings  using  a  range 
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Figure  6.  Wing  planform  view  and  view  from  wingtip  of  k  -  1  plane  of  C-H  grid. 


of  grid  sizes  from  coarse  to  very  fine.  He  found  that  for  a  medium  density 
grid,  surface  pressure  distributions  on  the  wing  surface  matched  well  with 
wind  tunnel  experiments.  For  grids  which  were  coarser  than  this,  however, 
he  found  that  the  expansion  peaks  in  the  surface  pressure  curves  were  more 
spread  out  and  the  absolute  value  of  their  minimum  pressure  coefficients 
were  reduced.  This  result  suggests  that  the  vortex  core  in  the  coarse  grid 
Euler  solution  was  larger  than  the  core  in  the  medium  grid  solution  and  the 
wind  tunnel  tests.  When  Rizzi  ran  simulations  on  the  very  fine  grid,  the 
absolute  value  of  the  minimum  pressure  coefficient  on  the  wing  was  much 
greater  than  for  the  wind  tunnel  tests,  and  the  region  of  the  wing  effected 
by  the  vortex  was  smaller.  This  result  implied  that  the  very  fine  grid 
allowed  a  vortex  which  was  smaller  than  in  the  physical  experiment. 

This  dependence  on  grid  density  for  sizing  of  the  vortex  core  is  a 
very  undesirable  feature  of  vortex  flowfield  analyses  using  Euler  solvers 
as  well  as  Navier-Stokes  codes.  Not  only  does  core  size  effect  the 
pressure  distribution  on  the  wing  surface,  but  if  the  cell  size  changes 
along  the  length  of  the  vortex,  the  resulting  change  in  core  size  can 
produce  an  axial  pressure  gradient.  This,  in  turn,  can  contribute  to  the 
development  of  an  axial  jet  or  velocity  defect,  depending  on  the  direction 
of  the  gradient.  Because  these  effects  are  numerical  rather  than  based  on 
the  physics  of  the  problem,  the  usefulness  of  Euler  solvers  for  vortex 
research  is  jeopardized. 

In  addition,  the  use  of  different  surface  boundary  conditions  is 
significant  for  problems  involving  the  formation  of  leading  edge  vortices 
on  wings  with  rounded  leading  edges .  Figure  7  shows  an  example  of  this  for 
a  wing  with  a  leading  edge  sweep  angle  of  71.2  degrees  and  a  NACA  0012 
airfoil.  The  figure  shows  a  crossflow  plane  cut  of  the  wing  with  crossflow 
velocity  vectors  drawn  for  two  flow  solutions  from  FL057  at  eighteen 
degrees  angle  of  attack.  The  solution  on  the  top  was  obtained  by 
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Crossflow  Velocity  Vectors 
Angle  of  Attack  is  18  Degrees 
p»urf«c  “  Pen  Boundary  Condition 


Crossflow  Velocity  Vectors 

Angle  of  Attack  is  18  Degrees 

Normal  Momemtum  Equation  Boundary  Condition 


Effect  of  surface  boundary  condition  on  vortex  rollup 


extrapolating  surface  pressures  from  the  boundary  cell  centers  to  the  wing 
surface  using  the  momentum  equation  normal  to  the  surface,  while  the  one  on 
the  bottom  came  from  setting  the  surface  pressure  equal  to  the  boundary 
cell  pressure.  As  can  be  seen,  the  normal  momentum  equation  boundary 
condition  prevented  the  vortex  from  developing.  These  numerical 
difficulties  with  the  formulation  of  Euler  solvers  and  thr  choice  of  grids 
limit  the  degree  to  which  they  can  be  used  in  vortical  flow  research. 
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CHAPTER  3 


MODEL  DEVELOPMENT 


3.1  VORTEX  CORE  MODEL 


A  model  for  the  vortex  core  is  derived  from  the  steady,  incompressible 
Navier-Stokes  equations  written  for  the  cylindrical  coordinate  system 
centered  on  a  vortex  core.  The  modeling  process  is  similar  to  that 
developed  by  Mager  [27]  and  Krause,  [28]  but  the  present  model  includes 
approximations  for  non-axisymmetric  wall  effects.  The  vortex  is  assumed  to 
be  slender  and  close  to  a  wing  surface  and  a  plane  of  symmetry.  The 
velocity  components  u,  v,  and  w  are  defined  along  the  x,  r,  and  & 
coordinates  which  correspond  to  axial,  radial,  and  circumferential 
directions.  Vith  these  assumptions,  the  Navier-Stokes  equations  can  be 
written: 
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(7a) 


dud v  v  l  dw 

—  +  —  +  -  + - 

dx  dr  r  r  dO 


du  du  uidu  Idp  +  fif  d2U  d2u  \  du  ^  132lt\ 

Udx  +  Vdr  +  rdd  pdx  p\dx2  dr 2  rdr  r2dd2J 


(7b) 


dv  dv  wdv  w2  1  p 

u - +  V - + - - - - — 

dx  dr  r  d6  r  p*r 


p(  d2v  d2v  1  dv  2  dw\ 

+  p\dx2  +  dr2* r dr  r 2  r2d62  r2 dd ) 


(7c) 


dW  dW  Wdw  VW  1  dp 

IL - +  U - + - +  -  - - — 

dx  dr  r  dO  r  pr  dO 

uf  d2w  d2w  \dw  w  1  d2w  2  dv\ 

where  p  is  the  pressure,  pis  the  density,  and  p  is  the  air  viscosity.  In 
a  manner  similar  to  boundary  layer  analysis,  a  small  parameter,  e  is 
defined  by  the  equation 


e  ■ 


L 


where  L  is  the  reference  axial  length  scale  of  the  vortex  and  <5r#/  is  the 
vortex  core  radius  at  x  -  L.  It  is  further  assumed  that  the  ratios  of 
radial  distances  and  velocities  to  their  axial  counterparts  are  also  of  the 
same  order  as  €  .  To  non-dimensionalize  equations  (7),  x  is  normalized  by 
L,  but  r  is  normalized  by  e  L.  In  this  way  the  very  small  radial 
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coordinate  is  scaled  to  be  of  the  sane  order  of  magnitude  as  the  axial 
coordinate.  In  a  similar  manner,  u  and  w  are  normalized  by  the  free  stream 

velocity,  U.  ,  but  v  is  normalized  by  eU.  .  Pressure  is  normalized  by  free 

2 

stream  dynamic  pressure,  pu.,  In  addition,  the  vortex  Reynolds  number  is 
defined  as 


Re 


PU  mL 
P 


When  the  non-dimensional  variables  are  substituted  into  equations  (7)  and 
common  multiples  of  L,  U.  and  e  are  divided  out  the  equations  can  be 
written 
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(8a) 
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(8c) 
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Next,  the  parameter  £  Is  assumed  to  scale  with  the  inverse  of  the 
square  root  of  the  Reynolds  number.  For  very  large  Reynolds  numbers,  It  is 
therefore  assumed  that  terms  of  order  £?  1/Re,  e,  and  l/(Re£)  are 
negligible.  After  eliminating  these  terms,  the  equations  of  motion 
become : 


da  dv  v  ^  1  dm 
dx  dr  r  € r  dO 


(9a) 


du  da  wdu  dp  1  fl  d  f  du\  1  d2u\ 

dx  dr  r  dQ  dx  Ree2\rdr\  dr)  r2d$2f 


(9b) 


W2  dp 

r  dr 


(9c) 


dw  dw  w dw  vw  l  dp 

U +  V  —  + - + - — 

dx  dr  er  dd  r  er  dQ 


+ 


1  d2U)\ 

r2de2) 


(9  d) 


The  above  equations  are  integrated  algebraically  by  defining  radial 
distributions  or  profiles  of  the  axial  and  circumferential  velocities 
within  the  core.  A  number  of  profiles  which  have  been  used  for  this 
purpose  are  described  in  Appendix  A.  Since  all  of  these  profiles  are 
axisymmetric,  derivatives  in  the  circumferential  direction  for  these 
profiles  are  zero,  and  the  equations  to  be  integrated  reduce  to 
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In  many  of  the  velocity  profiles,  such  as  the  Ranklne  vortex,  the  radial 
velocity,  v,  is  also  assumed  to  be  zero.  This  further  simplifies  the 
formulation.  For  the  purposes  of  the  present  work,  v  is  generally  not 
assumed  to  be  zero,  and  the  vortices  are  only  assumed  to  be  nearly 
axisymmetric ,  so  that  derivatives  in  the  circumferential  direction  are 
small  but  not  necessarily  zero. 

To  include  non- axisymmetric  wall  effects,  the  concept  of  image 
vortices  is  used.  The  non-axisymmetric  velocity  profiles  for  a  vortex  near 
a  wall  are  approximated  by  adding  to  the  axisymmetric  profiles  the 
velocities  induced  by  the  image  vortices.  In  order  to  facilitate 
evaluation  of  the  image  vortex  velocities,  a  wing-centered  Cartesian 
coordinate  system  is  defined  with  X,  Y,  and  Z  as  chordwise,  vertical,  and 
spanwise  coordinates  respectively.  This  Cartesian  coordinate  system  is  the 
same  as  that  used  by  the  Euler  solver.  Figure  8  illustrates  the 
relationship  between  the  vortex-centered  xrd  system  and  the  Cartesian 
system.  In  the  XYZ  system,  two  non-dimensional  distances  and  two  angles 
are  defined  as  follows: 
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Coordinate  system  relationships . 
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Figure  8  also  illustrates  these  distances  and  angles.  For  infinite  image 
vortices,  the  induced  v  and  w  velocities  in  the  vortex  coordinate  system 
are  given  by 
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See  Appendix  B  for  a  derivation  of  equations  (10). 

For  the  problems  considered  in  the  present  work,  tY  is  small  and  lz  is 
large  compared  with  1Y  Image  vortex  induced  u  velocities  are  therefore 
ignored.  The  image  vortex  induced  w  velocity  is  added  to  the  chosen 
axisymmetric  w  profile  to  obtain  the  total  w  velocity  distribution. 

With  the  circumferential  velocity  profile  specified,  equation  (9c)  is 
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integrated  from  the  vortex  center,  r  -  0,  to  the  edge  of  the  core,  r  -  6. 
After  removing  negligible  terms,  the  relationship  between  the  pressure  at 
the  core  edge  and  pressure  at  the  vortex  center  is  given  by 

r  2 

p(x  ,0)  =  p(x  ,6)  -  K  a  —  +  K  ,r2  C  t.  (11) 

a 

where  the  vortex  core  area,  a  and  the  circulation,  Tare  defined  as 

a  *  n6 2  and  r  =  2  n6wt 

respectively,  and  C*  is  the  wall  effect  parameter  given  by 


C* 


cos  <t>y  cos  <f>z 


COSQyCOStp  z 


(12) 


The  coefficients  K.  and  K±  depend  on  the  choice  of  the  circumferential 
velocity  profile.  In  addition,  for  cases  where  the  image  vortices  are  not 
infinite,  the  value  of  is  reduced  to  approximate  this  fact.  An  example 
of  the  evaluation  of  K,  and  Kx  for  the  case  of  a  Rankine  vortex  plus  the 
image  vortex  velocities  is  given  in  Appendix  D. 

The  v  velocity  profile  is  next  obtained  by  integrating  equation  (9a) 
in  the  radial  direction.  To  simplify  this  process,  the  equation  is  first 
integrated  using  the  axisymmetric  u  and  w  profiles  to  obtain  v  for  the 
axisymmetric  case.  The  non-axisymmetric  v  defined  by  equation  (10a)  is 
then  added  to  the  axisymmetric  v.  An  example  of  this  process  is  also  given 
in  Appendix  D.  With  the  profiles  for  the  three  components  of  velocity 
known,  equation  (9d)  is  integrated  over  the  area  of  the  vortex  core  to 
obtain: 
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(13) 


K  |  d  F  K  2du6  K  3 
r  dx  u6  dx  Ree2 au6 

where  the  values  of  Kj^  K2,  and  K3  depend  on  the  velocity  profiles  chosen. 
For  the  Rankine  vortex  circumferential  velocity  profile,  K2  and  K3  are  zero 
and  is  one.  In  this  case,  equation  (13)  gives 

a  «  r 


which  implies  that  core  area  increases  linearly  with  vortex  strength.  In 
the  velocity  profiles  defined  by  Mager,  K3  and  K2  are  not  constants  but 
depend  on  the  form  parameter  of  the  axial  velocity  profile.  This  form 
parameter  varies  upstream  of  bursting  from  1  for  a  uniform  axial  flow  to  0 
when  stagnation  occurs  at  the  vortex  center.  See  Appendix  D  for  an  example 
of  the  effect  this  has  on  equation  (13).  In  order  to  simplify  the  model, 
average  values  of  Kx  and  K2  are  used. 

For  such  complex  velocity  profiles  as  Mager' s,  it  is  more  difficult  to 
integrate  equation  (13).  Two  extreme  cases  are  considered,  each  of  which 
permits  (13)  to  be  integrated  easily.  First,  it  is  assumed  that  vortex 
strength  is  increasing  rapidly  and  the  viscous  term  of  (13)  can  be 
neglected.  For  this  case  the  equation  can  be  integrated  to  give 


where  K*  is  a  constant  of  integration.  On  the  other  hand,  consider  a  case 
where  vortex  strength  and  velocity  at  the  core  edge  are  constant.  For 
these  assumptions,  equation  (13)  would  be  solved  by 
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*a  [x-x0) 
Ree2  u6 


where  x0  is  a  constant  of  integration.  By  choosing  a  suitable  coefficient 
K5  to  replace  K3  the  above  expression  for  core  area  could  be  written 


a 


Based  on  these  two  results,  the  following  approximate  model  for  core  growth 
is  defined: 


a 


K. 


(14) 


where  Kg  is  a  core  growth  parameter  which  depends  on  the  velocity  profiles 
and  the  product  Re  e2. 

Next,  the  circumferential  velocity  profile  is  assumed  to  be  close  to 
that  of  a  Rankine  vortex,  so  the  radial  distribution  of  pressure  in  the 
vortex  core  can  be  approximated  by 


p(x,r)- p(x,0)  +  ^j  (p(x,<5)-p(x,0)).  (15) 

Using  equations  (11)  and  (15),  equation  (9b)  is  integrated  over  the  area  of 
the  core  to  give 
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where  the  subscript  0  designates  core  center  values.  It  is  worthwhile  to 
point  out  that  the  contribution  of  the  viscous  term  disappears  in  the 
integration  if  the  axial  velocity  profile  is  defined  as  in  Mager's  profile 
with  zero  slopes  at  the  center  and  the  edge  of  the  core. 

If  the  velocity  profiles  of  the  viscous  core  and  those  in  Euler 
solutions  are  similar,  then  the  major  discrepancy  of  the  Euler  formulation 
is  primarily  due  to  the  non-physical  estimation  of  core  size  and  core 
growth  rate.  Therefore,  the  non-physical  vortex  core  in  Euler  solutions 
can  be  replaced  with  a  core  that  is  derived  from  viscous  physics,  as  given 
in  equation  (16).  This  is  in  effect  equivalent  to  an  addition  of  the 
difference  between  the  model  and  the  Euler  solution  as  a  source  term  into 
the  Euler  formulation.  The  difference  can  be  given  as 
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where  the  subscript  E  denotes  flow  quantities  measured  in  the  Euler 
solution. 
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Equation  (17)  provides  momentum  source  terms  for  use  in  the  Euler 
solver.  Because  the  vortex  x  axis  is  not  alligned  with  the  wing  X  axis, 
components  of  the  source  term  for  each  cell  must  be  applied  to  the  X,  Y, 
and  Z  momencum  equations  in  the  Euler  solver. 


3.2  SHEAR  LAYER  MODEL 


The  viscous  shear  layer  next  to  a  highly  swept  wing  at  a  large  angle 
of  attack  is  too  complex  to  model  with  simple  physics.  The  flowfield 
involves  many  difficult  features  including  the  mechanism  for  leading  edge 
separation  and  creation  of  leading  edge  vortices,  reattachment  of  the  flow 
around  the  primary  vortex,  appearance  of  a  secondary  vortex  separation, 
etc.  As  a  simple  approximation  of  these  complexities,  an  analogy  is  made 
to  the  boundary  layer  on  a  flat  plate.  Boundary  layer  separation  can  be 
explained  as  being  caused  by  loss  of  total  pressure  in  the  boundary  layer 
due  to  viscosity.  This  same  concept  is  extended  to  explain  leading  edge 
separation  on  a  highly  swept  rounded  leading  edge.  Total  pressure  loss  in 
a  boundary  layer  is  approximated  as  being  of  such  a  magnitude  that  static 
pressure  at  the  boundary  layer  edge  is  equal  to  the  pressure  at  the  flat 
plate  surface.  For  the  rounded  leading  edge,  an  effective  shear  layer 
thickness  is  defined  such  that  total  pressure  losses  satisfy  this  same 
equality  of  static  pressures  at  the  shear  layer  edge  and  at  the  surface. 

The  magnitude  of  this  effective  thickness  will  be  different  than  for  a  flat 
plate  boundary  layer.  It  is  assumed,  however,  that  the  growth  profile  of 
the  effective  shear  layer  thickness  can  be  approximated  by  a  simple 
modification  to  the  well  known  equation  for  turbulent  boundary  layer  growth 
on  a  flat  plate  given  by  Schlichting  [38].  The  modified  equation  is 
written 
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Ktl 

(  fie  local)0'2 


(18) 


where  Relocal  is  the  local  Reynolds  number  based  on  the  distance  1  measured 
from  the  apex  of  the  delta  wing  and  <5S  is  the  effective  shear  layer 
thickness.  The  coefficient  is  an  effective  thickness  parameter  which 
equals  0.37  for  a  flat  plate,  but  which  is  much  greater  for  the  round 
leading  edge.  This  shear  layer  thickness  effect  is  incorporated  into  the 
Euler  solver  by  modifying  the  normal  momentum  equation  boundary  condition. 
Instead  of  extrapolating  pressure  from  the  cell  center  all  the  way  to  the 
wing  surface,  the  modified  boundary  condition  uses  the  relation 


wall 
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dp 
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where  pc<11  is  the  pressure  at  the  center  of  the  boundary  cell,  dp/dn  is  the 
pressure  gradient  in  the  direction  normal  to  the  surface,  and  Jn  is  the 
normal  distance  between  the  surface  and  the  cell  center. 

In  addition  to  corrections  for  the  surface  pressure,  the  analogy  with 
the  flat  plate  boundary  layer  is  used  to  calculate  momentum  source  terms 
for  the  boundary  cells.  These  source  terms  are  similar  to  those  calculated 
by  the  vortex  core  model.  They  represent  the  net  momentum  loss  experienced 
by  the  boundary  cells  due  to  the  shear  layer.  If  the  entire  shear  layer  is 
contained  within  the  boundary  cells,  the  skin  friction  coefficient,  Cf,  is 
a  measure  of  this  momentum  transfer  from  the  fluid  in  the  boundary  cell  to 
the  wing  surface.  Using  the  analogy  with  a  flat  plate  boundary  layer  once 
again,  the  equation  for  Cf  in  the  shear  layer  is  approximated  as 


Cf 


K, 


{  local 


1 0  •  2 


(20) 
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Kt  in  equation  (20)  is  a  wall  shear  stress  parameter  which  equals  0.074  for 
a  flat  plate  boundary  layer.  The  shear  stress  source  terms  are  added  into 
the  momentum  equations  in  a  direction  parallel  to  the  local  velocity  in 
each  boundary  cell. 


3.3  MODEL  IMPLEMENTATION 


The  vortex  core  model  is  incorporated  into  the  Euler  solver  as  a 
separate  subroutine  which  is  accessed  periodically.  The  control  parameters 
Kg,  K,,  and  Kt  are  input  into  the  program  by  the  user  and  used  unchanged 
throughout  the  solution  process.  The  exponents  Kx  and  K2  in  equation  (14) 
are  both  assumed  to  be  approximately  equal  to  1.  The  total  circulation  or 
vortex  strength  and  the  Euler  vortex  core  size  are  determined  from  the 
Euler  solution  at  axial  locations  along  the  vortex  corresponding  to 
convenient  locations  in  the  Euler  grid.  The  vortex  strength  is  obtained  by 
evaluating  the  vorticity  flux  through  the  surfaces  of  each  cell  and  summing 
them.  The  area  of  the  vortex  core  is  measured  from  the  cross  sectional 
areas  of  the  cells  which  contain  vorticity  of  strength  above  a  given 
threshold  level.  The  center  of  the  vortex  at  each  axial  station  is 
identified  as  the  cell  with  a  local  minimum  static  pressure  and  a  local 
maximum  vorticity  flux.  Locating  the  vortex  center  establishes  the  values 
of  C#  in  equation  (12).  The  axial  momentum  source  terms  at  each  axial 
position  along  the  vortex  are  then  obtained  from  equation  (17). 

When  the  momentum  source  terms  are  added  in,  the  equation  solved  by 
the  Euler  code  is  of  the  form 


where  the  S  term  represents  the  momentum  sources  from  the  core  model.  To 
avoid  destabilizing  the  solution,  source  terms  are  not  just  applied  to  the 
cell  at  the  vortex  center.  Terms  for  neighboring  cells  are  scaled  in 
accordance  with  their  distance  from  the  center  following  a  profile  similar 
to  the  axial  velocity  profiles  used  by  Mager.  In  order  to  reduce 
computational  costs,  these  corrections  are  updated  only  once  every  30  to  50 
iterations . 

The  new  boundary  condition  is  incorporated  into  the  Euler  solver  by 
modifying  the  surface  boundary  condition  subroutine.  The  desired  values  of 
effective  thickness  parameter  and  shear  stress  parameter  are  input  to  the 
program.  These  values  are  used  to  calculate  the  effective  viscous  layer 
thickness  for  each  boundary  cell  using  equation  (18).  The  pressure  at  the 
surface  is  then  defined  using  equation  (19) .  Momentum  source  terms  for  the 
boundary  cells  are  scaled  by  the  local  values  of  skin  friction  coefficient 
given  by  equation  (20). 
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CHAPTER  4 


CORE  MODEL  EVALUATION 


4.1  RECTANGULAR  GRID  EXPERIMENT 


A  simple  computational  experiment  was  performed  which  simulated  a 
vortex  near  an  infinite  flat  plate  and  a  plane  of  symmetry.  The  experiment 
had  three  purposes.  The  first  of  these  was  to  identify  some  of  the 
characteristics  of  Euler  codes  which  cause  problems  in  vortex  simulations. 
The  second  was  to  evaluate  some  of  the  assumptions  and  approximations  used 
in  developing  the  vortex  core  model.  The  third  purpose  was  to  provide  a 
simple  computational  problem  for  the  first  tests  of  the  vortex  core  model. 

A  computational  grid  was  defined  in  the  shape  of  a  rectangular  box. 

The  flowfield  within  the  box  was  initialized  with  a  uniform  flow  parallel 
to  the  top,  bottom,  and  sides  of  the  box.  To  this  was  added  a  vortex  whose 
axis  ran  parallel  to  the  uniform  flow.  This  vortex  flowfield  was  also  used 
as  an  inflow  boundary  condition  on  the  upstream  end  of  the  box.  A  solid 
surface  boundary  condition  was  imposed  on  the  bottom  surface  of  the  box  and 
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one  of  the  sides  of  the  box  was  made  a  plane  of  symmetry.  The  other  three 
box  faces  were  given  farfield  boundary  conditions. 

The  motivation  of  this  study  was  to  isolate  the  difficult 
computational  issues  arising  from  the  complicated  geometries  and  flowfields 
associated  with  delta  wings.  Such  factors  as  leading  edge  separation, 
vortex  rollup,  and  the  pressure  fields  associated  with  wing  airfoil  shapes 
were  completely  excluded  from  the  test.  Other  factors  such  as  vortex 
strength,  location,  and  core  size  could  be  controlled  individually.  The 
two  different  Euler  solvers  discussed  in  Chapter  2  were  used  in  the  tests. 
Both  gave  essentially  the  same  results. 

The  first  test  involved  initializing  the  flow  field  with  a  potential 
vortex.  In  order  to  avoid  a  singularity  at  the  vortex  center,  the  vortex 
was  centered  on  a  grid  node  so  that  induced  velocities  at  the  cell  centers 
were  finite  for  FL057  and  vice  versa  for  ARC3D.  This  caused  the  entire 
flow  field  to  be  initially  irrotational .  This  velocity  distribution  could 
also  be  regarded  as  that  from  a  Rankine  vortex  with  a  core  diameter  equal 
to  the  width  of  a  single  grid  cell.  The  velocity  fields  induced  by  three 
image  vorticies  were  also  added  to  the  initialization  so  that  normal 
velocities  at  the  flat  plate  and  the  plane  of  symmetry  were  zero.  When  the 
Euler  code  was  run,  the  action  of  the  numerical  viscosity  in  the  flow 
solver  caused  the  vortex  core  to  enlarge  to  cover  a  region  eight  to  ten 
cells  in  diameter.  This  growth  in  core  size  occurred  within  the  first  2 
cells  downstream  of  the  upstream  boundary.  Total  swirl  or  vortex  strength 
was  conserved,  so  maximum  circumferential  velocity  in  the  vortex  decreased 
dramatically  in  this  relatively  short  downstream  distance.  The  reduction 
in  maximum  circumferential  velocity  was  accompanied  by  an  increase  in 
pressure  along  the  central  axis  of  the  vortex.  Since  the  upstream  boundary 
condition  was  inflow,  the  velocity  profile  there  did  not  change  with 
iterations.  With  a  potential  vortex  at  the  upstream  boundary  and  a  more 
diffuse  core  just  a  few  cells  downstream,  a  strong  pressure  gradient  was 
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created  which  was  opposite  to  the  direction  of  the  axial  flow  in  the  vortex 
core.  This  adverse  pressure  gradient  caused  a  deficit  in  the  axial  flow 
velocities  in  the  core  region  and,  in  fact,  the  flow  slowed  to  near 
stagnation  in  the  cells  just  downstream  of  the  upstream  boundary. 

To  avoid  the  near-bursting  condition  caused  by  the  potential  vortex 
initial  condition,  a  more  "natural"  circumferential  velocity  profile  was 
used.  This  velocity  distribution  was  obtained  from  the  original  potential 
vortex  solution  at  a  location  5  to  8  cells  downstream  of  the  point  where 
the  core  grows  to  its  equilibrium  size.  The  shape  of  this  profile  was 
determined  by  the  numerical  viscosity,  but  it  was  very  similar  to  Mager's 
analytic  velocity  profile  (see  Appendix  A).  Figure  9  compares  radial 
distributions  of  circumferential  velocity  and  circulation  for  Mager's 
profile  and  for  a  Rankine  vortex  with  a  typical  profile  from  the  Euler 
solution.  Note  that  the  Euler  results  are  extremely  close  to  the  Mager 
profile . 

The  solutions  obtained  with  this  Mager- like  upstream  boundary 
condition  are  quite  well  behaved.  Such  flow  quantities  as  core  diameter, 
total  circulation,  and  axial  velocity  are  essentially  constant  down  the 
length  of  the  vortex.  The  only  quantity  which  varies  significantly  in  the 
axial  direction  is  the  distance  of  the  vortex  from  the  flat  plate  and  the 
symmetry  plane.  The  image  vortex- induced  velocity  resulting  from  the 
proximity  of  the  vortex  to  a  solid  wall  causes  the  distance  of  the  vortex 
from  the  plane  of  symmetry  to  increase  with  increasing  downstream  distance. 
The  wall  effect  also  causes  the  vortex  to  move  slightly  closer  to  the  flat 
plate  as  it  moves  downstream.  The  magnitudes  of  the  self- induced 
velocities  in  the  Euler  solution  agree  very  well  with  what  is  predicted 
using  the  concept  of  image  vortices.  A  very  slight  increase  in  minimum 
core  pressure  was  noted  as  the  vortex  moved  closer  to  the  flat  plate. 

Figures  10,  11,  12,  and  13  show  crossflow  plane  streamlines,  velocity 
vectors,  vorticity  contours,  and  pressure  contours  for  a  typical  flowfield 
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Figure  10.  Rectangular  grid  study  crossflow  streamlines.  Flat  plate  surface 
is  at  the  bottom  and  plane  of  symmetry  is  at  right. 
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Figure  11.  Rectangular  grid  study  crossflow  velocity  vectors.  Flat  plate 
surface  is  at  the  bottom  and  plane  of  symmetry  is  at  right. 


Figure  13.  Rectangular  grid  study  < 
plate  surface  is  at  the  bottom  and  ] 
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at  a  point  halfway  down  the  length  of  the  vortex.  Tick  marks  on  the  frame 
of  the  two  latter  plots  correspond  to  individual  cells  on  the  uniform 
rectangular  grid.  Note  that  the  region  of  reduced  static  pressure  spans 
approximately  ten  cells.  The  region  of  significant  vorticity  is  also  about 
ten  cells  in  diameter.  Outside  this  rotational  core,  however,  the  rest  of 
the  flowfield  is  very  nearly  irrotational .  This  result  suggests  that  the 
use  of  irrotational  image  vortices  in  the  vortex  core  model  is  acceptable. 
As  a  further  evaluation  of  the  use  of  irrotational  image  vortices  in  the 
core  model.  Figure  14  shows  a  crossflow  streamline  plot  for  a  Hager  profile 
with  image  vortex  velocities  added  to  it.  The  value  of  the  stream  function 
used  for  this  plot  is  adjusted  so  that  the  streamline  labeled  0  is  the  core 
edge.  Note  the  very  good  agreement  between  this  plot  and  Figure  10. 

The  relationship  in  the  core  model  for  the  variation  of  core  pressure 
with  distance  from  walls  can  be  evaluated  by  comparing  it  to  pressure 
variations  in  the  Euler  solutions.  Figure  15  plots  vortex  core  pressure 
for  two  Euler  solutions  as  a  function  of  the  wall  proximity  parameter  C#. 
The  solid  line  on  the  plot  represents  the  pressure  variation  predicted 
using  the  vortex  core  model  equation  (11).  Equation  (11)  appears  to  be  a 
very  good  approximation  for  the  effect  of  walls  on  core  pressure. 

Several  tests  were  made  to  evaluate  grid  effects  in  the  flow 
solutions.  The  first  test  was  a  study  of  the  effect  of  the  overall  size  of 
the  grid  on  the  flow  solution.  It  was  observed  early  in  the  research  that 
if  the  grid  was  so  small  that  the  image  vortex  velocity  caused  the  vortex 
to  move  into  the  cells  adjacent  to  the  lateral  far  field  boundary,  the 
boundary  conditions  used  could  not  keep  the  vortex  core  intact.  For  a  grid 
with  18  cells  in  the  X  direction,  a  Y-Z  plane  which  was  32  cells  high  and 
48  cells  wide  seemed  to  give  good  results,  with  the  vortex  never  getting 
clos  r  to  the  far  field  than  Y  -  11  and  Z  -  23.  When  this  grid  was  reduced 
in  size  to  24  by  36,  a  0.3%  increase  could  be  seen  in  the  minimum  pressures 
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Figure  15.  Variation  of  minimum  core  pressure  with  wall  proximity  factor,  C( 


along  the  vortex  core  and  the  image  vortex- induced  velocity  of  the  vortex 
core  was  10%  less.  However,  when  the  grid  size  was  increased  to  40  x  60, 
the  difference  in  minimum  core  pressures  was  less  than  0.1%  and  the  lateral 
motion  of  the  core  showed  no  difference  from  that  on  the  32  x  48  grid.  A 
study  was  also  made  of  the  effect  of  stretching  the  outer  cells  to  place 
the  far  field  boundaries  even  further  away.  For  the  18  x  32  x  48  grid, 
stretching  of  the  outer  6  cells  in  both  the  Y  and  Z  directions  had  little 
effect  on  the  solutions. 

In  order  to  evaluate  the  effects  of  cell  aspect  ratio  on  the  quality 
of  the  flow  solution,  several  test  runs  were  made  with  the  cells  elongated 
or  shortened  by  changing  the  grid  spacing  in  the  X  or  streamwise  direction. 
In  general,  the  effect  of  stretching  or  shortening  the  cells  was  negligible 
over  a  grid  spacing  from  .5  to  2.  in  the  X  direction  with  the  Y  and  Z  grid 
spacing  constant  at  1. 

The  effect  of  grid  skewing  on  the  quality  of  the  solution  was  also 
investigated.  Since  the  image  vortex  velocities  caused  the  vortex  to  move 
at  an  angle  from  the  axial  direction,  the  grid  was  skewed  at  this  angle  to 
keep  the  vortex  in  the  cell  with  the  same  k  index  in  each  successive  i 
plane.  The  cells  in  the  region  between  the  vortex  and  the  plane  of 
symmetry  were  stretched  as  necessary  to  keep  the  plane  of  symmetry  parallel 
to  the  free  stream  flow  direction.  Figure  16  shows  three  views  of  one  of 
these  skewed  grids,  the  grid  labeled  Grid  A.  The  solutions  obtained  on 
this  grid  for  the  same  number  of  iterations  differed  less  than  0.04%  from 
those  obtained  on  the  uniformly  spaced  grid,  an  almost  indistinguishable 
difference.  However,  the  calculation  converged  faster  on  the  skewed  grid, 
so  the  slight  error  which  appears  to  exist  between  the  two  may  be  mostly  a 
measure  of  convergence.  This  result  suggests  that  a  reasonable  amount  of 
grid  skewness  is  tolerable  in  flow  calculations  around  more  complex 
geometries . 
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Figure  16.  Skewed  and  stretched  grids. 
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Grid  A  in  Figure  16  has  a  much  finer  grid  spacing  in  the  region  where 
the  vortex  core  was  placed,  but  it  is  stretched  in  the  outer  regions  so  the 
far  field  is  kept  at  the  same  distance  as  it  was  on  the  uniform  grid.  As 
expected,  a  finer  grid  allowed  a  smaller  vortex  core  since  the  core  spread 
over  the  same  number  of  cells  as  in  the  uniform  grid  case.  In  order  to 
demonstrate  the  error  this  cell  size  effect  could  cause  in  a  flow  solution, 
a  grid  was  generated  which  had  a  grid  spacing  of  one  half  unit  in  the  area 
of  the  vortex  at  the  upstream  end,  but  a  spacing  of  2  units  at  the 
downstream  end  of  the  domain.  Three  views  of  the  grid  are  also  shown  in 
Figure  16  and  labeled  Grid  B.  The  grid  spacing  in  the  streamwise  or  X 
direction  was  one  unit  throughout  the  domain,  so  the  cells  expanded  from  .5 
x  .5  to  2  x  2  in  a  distance  of  18  units.  This  rate  of  grid  expansion  is 
much  less  than  that  of  many  grids  commonly  used  for  computations  of  flows 
around  wings  (see  Rizzi  [19]  for  example).  However,  when  a  vortex  with  a 
core  of  the  minimum  size  resolvable  by  the  Euler  code  on  a  given  region  of 
the  grid  was  caused  to  flow  into  a  part  of  the  grid  where  the  cell  size  was 
increasing,  the  vortex  core  also  enlarged.  Since  the  flow  solver  preserved 
the  total  swirl  of  the  vortex,  this  non-physical  growth  of  the  core  caused 
reduced  circumferential  velocities  and  increased  pressures  in  proportion  to 
the  enlargement  of  the  cells.  Figure  17  compares  side  views  of  the  two 
grids  with  velocity  vector  arrows  drawn  for  their  respective  flow 
solutions.  In  the  flowfield  for  Grid  B,  it  is  apparent  that  the  adverse 
pressure  gradient  within  the  vortex  core  has  created  a  reversal  in  the 
axial  flow  within  the  core  similar  to  that  which  characterizes  vortex 
bursting.  The  test  case  run  with  identical  flow  conditions  and  boundary 
conditions  on  Grid  A  shows  only  a  slight  reduction  in  axial  flow 
velocities.  Figure  18  plots  the  axial  variation  of  pressure  for  the  two 
solutions  of  Figure  17. 
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Figure  17.  Side  views  of  velocity  vectors  in  the  vertical  plane  containing 
the  vortex  axis  for  two  flow  solutions  with  identical  initial  and  boundary 
conditions  but  different  grids. 


Two  general  conclusions  can  be  drawn  from  the  results  of  these  tests. 
First,  the  adverse  effects  of  grid  density  and  numerical  viscosity  on  flow 
solutions  containing  vortices  were  shown  to  be  potentially  very  serious,  to 
the  point  of  causing  bursting.  Second,  many  of  the  assumptions  used  in 
deriving  the  vortex  core  model  were  shown  to  be  reasonable  and  compatible 
with  the  Euler  formulation. 

4.2  BURSTING  CONTROL  DEMONSTRATION 

The  vortex  core  model  should  be  able  to  correct  for  or  reduce  the 
effects  of  grids  and  numerical  viscosity  in  two  ways.  If  axial  velocity 
defects  or  bursting  are  caused  by  an  expanding  grid,  the  vortex  core  model 
should  be  able  to  correct  for  that  effect  when  the  physical  vortex  core  is 
modeled  with  no  rate  of  growth.  On  the  other  hand,  if  the  grid  is  uniform 
and  little  or  no  axial  velocity  defect  is  present  in  the  flow  solution, 
then  by  modeling  a  physical  core  with  a  rapid  rate  of  growth,  the  vortex 
core  model  should  be  able  to  produce  bursting. 

Both  these  capabilities  were  successfully  demonstrated.  Figure  19 
shows  side  views  of  velocity  vectors  in  the  vertical  plane  containing  the 
vortex  for  two  different  flow  solutions  on  Grid  A.  The  solution  on  the 
left  was  obtained  using  an  Euler  solver  alone.  The  one  on  the  right  came 
from  an  Euler  solver  with  the  vortex  core  model.  For  this  test,  the 

physical  core  was  modeled  as  growing  very  rapidly  by  giving  the  vortex 

growth  parameter  a  large  value.  The  vortex  bursting  produced  by  the  model 
is  apparent. 

Figure  20  shows  a  demonstration  of  the  model's  ability  to  eliminate 
grid- induced  vortex  bursting  on  Grid  B.  The  flow  solution  with 

grid- induced  bursting  on  the  left  is  for  an  Euler  solver  alone,  while  the 


solution  on  the  right  is  for  the  Euler  solver  with  core  model.  For  this 
test,  the  physical  core  size  was  arbitrarily  set  to  the  same  constant  value 
down  the  length  of  the  vortex.  Since  no  axial  pressure  gradient  should 
develop  if  core  area  is  constant,  the  model  generated  source  terms  which 
eliminated  the  grid- induced  bursting. 
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CHAPTER  5 


RESULTS  AND  DISCUSSION 


5 . 1  DELTA  WING  BURSTING  STUDY 


The  vortex  core  model  was  tested  to  assess  its  ability  to  predict  the 
onset  of  bursting  on  four  different  delta  wings.  Grids  with  H-0  topology 
were  used  around  flat  plate  delta  wings  with  sweep  angles  varying  from  55 
to  76  degrees.  All  grids  used  were  57  x  21  x  33  with  25  x  33  points  on  the 
wing  surface.  Figure  21  shows  planform  views  of  the  right  half  of  each 
wing  used  for  the  test.  Mach  number  for  all  these  tests  was  0.2. 

The  first  wing  tested  was  swept  65  degrees.  Flow  solutions  were  first 
obtained  for  this  wing  using  the  Euler  solver  alone.  It  was  discovered 
that  at  angles  of  attack  above  25  degrees,  some  bursting  was  present  in  the 
Euler  solutions.  This  was  not  surprising,  since  cell  sizes  in  the  grid 
increased  in  the  streamwise  direction,  forcing  the  vortex  core  in  the  Euler 
solution  to  grow.  The  locations  of  the  bursting  points  did  not,  however 
agree  with  available  wind  tunnel  data  [4]  for  this  wing.  This  suggests 
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vl«wa  of  the  right  halves  of  four  delta  wings  used  In  the 


that  the  rate  of  core  growth  induced  by  the  grid  was  not  the  same  as  that 
of  the  physical  vortex  in  the  wind  tunnel  tests. 

For  the  core  model  tests,  a  value  for  Kg  was  chosen  which  gave  the 
same  core  area  at  the  wing  trailing  edge  as  that  in  the  Euler  solution. 

For  this  situation,  the  difference  between  the  Euler  core  and  the  model 
core  was  limited  to  differences  in  the  core  growth  profile,  but  not  in  the 
final  core  size.  Appropriate  values  of  Ka  and  were  determined  by  trial 
and  error.  The  values  calculated  for  the  Rankine  vortex  and  Mager's 
profiles  were  used  as  initial  guesses,  but  these  proved  to  be  too  large. 
When  these  large  values  were  used,  the  solution  was  destabilized  and  caused 
to  diverge  rapidly.  As  smaller  values  of  Ka  and  Kt  were  tried,  a  range  was 
found  where  they  were  not  so  small  that  they  did  not  produce  an  effect,  but 
not  so  large  that  they  caused  instability.  Well  within  the  bounds  of  this 
range,  values  were  found  which  allowed  the  model  to  generate  source  terms 
which  brought  the  bursting  locations  in  the  resulting  solution  into  very 
close  agreement  with  the  wind  tunnel  data.  Figure  22  shows  a  comparison 
between  bursting  locations  for  the  Euler  alone  solutions  and  the  results 
from  the  Euler  code  with  core  model.  Wind  tunnel  data  for  the  same  wing 
are  also  plotted.  As  can  be  seen,  the  agreement  between  the  wind  tunnel 
data  and  the  model  results  is  very  good  over  the  entire  range  of  angles  of 
attack,  except  at  the  very  highest  angles.  Here  the  Euler  solver  with 
model  could  never  quite  get  bursting  to  progress  all  the  way  to  the  apex. 
This  trend  is  visible  in  the  Euler  alone  solution  as  well,  and  was  common 
in  all  the  other  delta  wing  bursting  tests  made. 

At  the  lower  angles  of  attack,  the  method  experienced  some 
instability.  This  was  overcome  by  increasing  the  number  of  iterations  of 
the  Euler  solver  between  updates  of  the  core  model  momentum  source  terms. 
The  source  terms  were  applied  at  each  iteration,  but  their  values  were 
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updated  only  periodically.  For  angles  of  attack  above  25  degrees,  source 
terms  were  updated  every  50  iterations  of  the  Euler  solver.  For  20  and  23 
degrees  angle  of  attack,  updates  were  done  only  once  every  150  iterations. 
In  both  cases,  as  the  flow  solutions  progressed,  the  magnitude  of  the 
source  terms  decreased  with  each  update.  In  the  final  converged  solutions, 
the  momentum  sources  were  usually  about  10%  or  less  of  their  initial 
values . 

The  next  test  of  the  core  model  was  to  determine  if  the  same  values 
for  Kg,  Ka,  and  would  produce  equally  good  results  on  delta  wings  with 
much  different  planforms.  A  75  degree  delta  wing  was  first  tested.  Figure 
23  shows  bursting  location  curves  for  this  wing  for  Euler,  Euler  with  core 
model,  and  wind  tunnel  results.  Agreement  with  the  wind  tunnel  data  of 
Wentz  [4]  is  better  than  Euler  alone  results,  but  not  as  good  as  for  the  65 
degree  delta  wing.  When  the  wind  tunnel  data  of  Payne  and  Nelson  for  this 
wing  is  plotted,  however,  the  agreement  is  improved. 

Tests  were  also  made  on  a  55  degree  delta.  Figure  24  plots  these 
results.  Once  again,  the  agreement  between  Euler  with  model  results  and 
wind  tunnel  data  is  better  than  for  the  Euler  alone  solutions,  though  not 
as  good  as  for  the  65  degree  wing. 

In  order  to  compare  results  from  the  core  model  with  a  Navier- Stokes 
solution,  a  fourth  delta  wing  was  tested.  This  wing  had  an  aspect  ratio  of 
1  and  a  sweep  angle  of  approximately  76  degrees.  Figure  25  compares  the 
resulting  bursting  location  versus  angle  of  attack  curves  with  the  single 
Navier-Stokes  result  [39],  Of  significance  in  evaluating  the  usefulness  of 
the  present  method  is  the  fact  that  all  of  the  Euler  and  Euler  with  model 
solutions  obtained  for  this  test  required  less  total  computation  time  than 
the  single  Navier-Stokes  result. 

Figure  26  shows  spanwise  surface  pressure  coefficient  distributions  at 
the  midpoint  of  the  root  chord  of  the  aspect  ratio  1  delta  wing  at  20.5 
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Figure  26.  Spanvise  pressure  coefficient  distributions  at  X/C  -  0.5  for  an 
AR  -  1  delta  wing  at  20.5  degrees  angle  of  attack. 
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degrees  angle  of  attack.  At  this  angle  of  attack,  the  vortex  does  not 
burst  over  any  part  of  the  wing.  The  Cp  curve  predicted  by  the  Euler 
solver  with  the  core  model  is  compared  with  curves  from  the  Euler  solver 
alone  and  from  Navier-Stokes  simulations  [39]  and  wind  tunnel  tests  [2], 
Note  that  neither  of  the  Euler  solvers  captured  the  secondary  vortex 
influence  on  the  Cp  curve  which  is  evident  in  the  wind  tunnel  data  and 
Navier-Stokes  solutions.  However,  the  vortex  core  model  does  allow  the 
modified  Euler  code  to  predict  the  peak  of  the  Cp  curve  in  closer  agreement 
with  experiment.  Although  bursting  does  not  occur  in  this  flow  solution, 
the  core  model  still  has  an  effect  on  the  axial  velocities  in  the  core.  As 
axial  velocities  are  decreased,  the  core  expands  slightly,  reducing 
circumferential  velocities  and  increasing  core  pressure.  Figure  27 
compares  velocity  vectors  in  a  nearly  horizontal  plane  which  is  parallel  to 
the  Z  axis  and  contains  the  vortex  axis  for  the  same  wing  at  40  degrees 
angle  of  attack.  The  result  from  the  Euler  code  alone  shows  the  bursting 
point  at  an  axial  station  about  70%  of  the  root  chord  downstream  of  the 
wing  apex.  The  Euler  code  with  core  model  gives  the  bursting  point  at 
about  30%  root  chord,  which  does  agree  better  with  Navier-Stokes  results 
[ 39  j .  Figure  28  compares  the  corresponding  Cp  curves  for  the  two  solutions 
at  50%  root  chord.  The  vortex  at  that  station  has  burst  in  the  Euler  with 
model  solution,  but  not  in  the  Euler  alone  solution.  Note  the  flattening 
of  the  Cp  curve  when  bursting  has  occurred. 


5.2  BOUNDARY  CONDITION  COMPARISONS 


The  shear  layer  model  surface  boundary  condition  was  tested  on  the 
arrow  wing  shown  in  Figure  6.  As  a  basis  for  comparison,  test  runs  were 
first  made  using  FL057  with  two  choices  of  surface  boundary  conditions. 
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The  first  of  these  extrapolates  surface  pressure  from  the  boundary  cell 
center  using  the  momentum  equation  normal  to  the  wing  surface.  The  second 
simply  equates  pressure  at  the  surface  to  pressure  at  the  cell  center. 
Surface  pressure  plots  for  a  range  of  angles  of  attack  for  the  two  boundary 
conditions  are  compared  with  wind  tunnel  data  [40]  in  Figure  29.  The 
tendency  for  the  momentum  boundary  condition  to  delay  leading  edge 
separation  is  apparent.  Likewise,  the  vortex  appears  to  develop  too 
quickly  when  using  the  simple  boundary  condition  which  equates  surface 
pressure  to  boundary  cell  pressure.  Figure  30  shows  results  using  the  new 
boundary  condition  with  the  shear  layer  model  at  four  degrees  angle  of 
attack  for  a  range  of  values  of  the  effective  shear  layer  thickness 
parameter,  Kt.  Wind  tunnel  results  are  shown  for  comparison.  Reynolds 
number  is  6  million  for  both  the  computer  model  and  the  wind  tunnel  tests. 

A  value  of  2.1  for  Kt  appears  to  give  the  best  results. 

Figure  31  compares  surface  pressure  plots  from  wind  tunnel  tests  with 
those  for  the  shear  layer  model  surface  boundary  condition  for  a  range  of 
angles  of  attack.  Kt  was  held  constant  at  2 . 1  for  this  simulation.  Note 
that  the  new  boundary  condition  results  match  the  wind  tunnel  data  better 
than  results  for  either  of  the  old  boundary  conditions,  especially  for 
angles  of  attack  of  eight  degrees  and  less.  At  10  degrees,  a  difference 
can  be  seen  between  the  location  of  the  Euler  with  model  vortex  and  the 
wind  tunnel  vortex.  This  difference  may  be  due  to  the  crude  approximation 
of  the  viscous  effects  in  the  shear  layer  used  in  developing  the  present 
model . 

For  the  preceding  tests,  the  shear  stress  parameter,  Ks  was  left  equal 
to  zero.  For  another  set  of  tests,  Ks  was  set  equal  to  0.074  and  used  in 
the  Euler  solver  with  the  shear  layer  model  on  two  different  cropped  delta 
wings  at  relatively  large  angles  of  attack.  Both  wings  had  sweep  angles  of 
65  degrees  and  taper  ratios  of  0.15.  The  first  wing  used  a  NACA  0012 
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airfoil.  When  a  flow  solution  was  obtained  for  this  wing  at  30  degrees 
angle  of  attack,  a  secondary  vortex  was  seen  in  the  flowfield  between  the 
primary  vortex  and  the  leading  edge.  Figure  32  shows  crossflow  velocity 
vectors  at  X/Croot  -  0.5  for  this  result  compared  with  results  from  an  Euler 
solver  alone  for  the  same  wing  and  flow  conditions.  The  secondary  vortex 
is  made  especially  visible  by  the  contrast  with  the  Euler  alone  solution. 

In  the  second  test,  a  delta  wing  with  the  same  planform  but  a  NACA 
65A005  airfoil  was  tested  at  20  degrees  angle  of  attack.  Ks  was  again  set 
to  0.074.  Figure  33  shows  surface  pressure  contours  for  this  result 
compared  Euler  alone  results  for  the  same  wing  and  flow  conditions.  Note 
that  the  secondary  vortex  causes  the  primary  vortex  to  move  inboard  and 
that  a  small  local  trough  of  pressure  coefficient  develops  under  the 
secondary  vortex.  Figure  34  shows  crossflow  velocity  vectors  for  these 
conditions.  The  secondary  vortex  is  just  barely  visible  as  only  a  single 
row  of  three  arrows  a  little  inboard  of  the  leading  edge.  Figure  35  shows 
velocity  vectors  on  the  surface  of  the  wing  for  the  Euler  alone  solution 
and  the  Euler  with  model  results.  The  secondary  separation  line  is 
visible,  as  is  the  area  of  inboard  flow  under  the  secondary  vortex. 
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Figure  33.  Surface  pressure  contours  for  a  65  degree  delta  wing  with  a  NACA 
65A005  airfoil  at  20  degrees  angle  of  attack. 


Figure  34.  Crossflow  velocity  vectors  at  X/Croot  -  0.87  for  a  65  degree  delta 
wing  with  a  NACA  65A005  airfoil  at  20  degrees  angle  of  attack. 


Planform  view  of  surface  velocity  vectors  on  a  65 
NACA  65A005  airfoil  at  20  degrees  angle  of  attack 


CHAPTER  6 


CONCLUSIONS 


A  vortex  core  model  has  been  derived  which  introduces  the  effects  of 
physical  viscosity  for  vortex  cores  into  computer  solutions  of  the 
3 -dimensional  Euler  equations.  The  formulation  of  this  model  allows  it  to 
replace  the  effects  of  numerical  viscosity  in  the  solution  with  those  based 
on  physical  viscosity.  The  model  has  been  used  to  predict  the  onset  and 
location  of  vortex  bursting  on  a  variety  of  delta  wings.  In  each  case,  the 
Euler  code  with  the  core  model  predicted  bursting  locations  as  a  function 
of  angle  of  attack  which  were  in  much  better  agreement  with  wind  tunnel 
data  and  Navier-Stokes  simulations  than  the  results  for  an  Euler  code 
alone.  The  same  values  of  the  model  control  parameters  produced  reasonably 
good  results  for  all  the  tests.  Such  good  performance  for  a  variety  of 
wings  and  angles  of  attack  without  any  adjustments  of  the  model  control 
parameters  is  a  strong  argument  that  the  approximations  made  in  deriving 
the  model  are  acceptable,  and  that  the  model  is  useful.  The  core  model 
also  produced  favorable  effects  on  wing  surface  pressure  distributions 
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which  agree  with  wind  tunnel  results  and  Navier-Stokes  simulations  better 
than  results  from  the  Euler  code  alone.  This  improvement  is  seen  whether 
or  not  bursting  is  present  in  the  solution. 

Another  model  which  describes  viscosity  effects  in  the  viscous  shear 
layer  near  highly  swept  wings  has  also  been  derived.  The  model  uses  an 
analogy  to  the  boundary  layer  on  a  flat  plate  to  make  a  correction  to  the 
pressure  on  the  wing  surface.  This  pressure  correction  accounts  for  total 
pressure  losses  due  to  viscous  dissipation  in  the  shear  layer,  and 
contributes  to  the  leading  edge  separation  which  rolls  up  into  a  leading 
edge  vortex.  The  model  also  calculates  momentum  source  terms  for  the  grid 
cells  next  to  the  wing  surface  to  account  for  the  net  momentum  loss  in 
these  cells  due  to  viscosity.  The  new  boundary  condition  has  been  tested 
on  several  highly  swept  wings  with  rounded  leading  edges  for  a  range  of 
angles  of  attack.  Surface  pressure  distributions  for  these  tests  were 
found  to  be  in  better  agreement  with  wind  tunnel  tests  of  the  same  wing 
geometry  than  results  for  an  Euler  code  with  conventional  boundary 
conditions.  At  higher  angles  of  attack,  a  secondary  separation  and 
secondary  vortex  could  be  seen  in  the  flow  solution.  This  feature  makes 
these  flow  solutions  more  similar  to  Navier-Stokes  solutions  and  wind 
tunnel  tests  than  any  solutions  from  Euler  solvers  alone. 

The  two  models  for  vortex-related  viscosity  effects  which  have  been 
developed  in  the  present  work  can  be  extended  in  many  ways  to  improve  their 
capabilities.  The  success  of  these  models  and  the  great  need  for  the 
capability  they  can  provide  make  such  extensions  very  attractive. 

The  next  logical  step  for  the  vortex  core  model  is  to  extend  it  to 
include  compressibility  effects  so  the  model  could  be  used  for  higher  Mach 
number  problems.  An  extension  of  the  model  to  include  unsteady  flow  could 
also  be  done.  This  might  have  applications  for  time-accurate  simulations 
of  vortex  bursting  used  to  determine  aircraft  stability  derivatives. 

Higher  order  approximations  for  the  vortex  core  growth  model  could  also  be 
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developed.  This  could  even  be  formulated  as  a  differential  equation  which 
is  integrated  numerically  whenever  the  model  sources  terms  are  updated. 

The  shear  layer  model  could  also  be  extended  to  include 
compressibility  and  unsteady  effects.  The  first  step,  however,  should  be 
to  use  a  model  for  the  shear  layer  thickness  and  shear  stress  which 
accounts  for  geometry  effects  such  as  leading  edge  curvature  and  sweep 
angle  and  the  presence  of  a  flow  reattachment  line.  This  should  permit  the 
model  to  be  used  for  a  greater  variety  of  wing  shapes  without  changing  the 
control  parameters . 
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APPENDIX  A 


VORTEX  VELOCITY  PROFILES 


A  great  number  of  vortex  velocity  profiles  have  been  defined  as  solutions 
to  simplified  versions  of  the  Navier- Stokes  equations  or  as  approximations  to 
the  velocity  distributions  measured  in  real  vortices.  The  most  important  of 
these  will  be  described  in  a  generally  chronological  order. 

The  point  vortex  is  a  two-dimensional  circulating  flow  which  solves 
Laplace's  equation.  The  line  vortex  is  a  three-dimensional  extension  of  the 
point  vortex.  Both  are  sometimes  referred  to  as  potential  vortices  because 
their  flowfields  are  irrotational  everywhere  but  at  their  exact  center.  The 
line  vortex  velocity  profile  is  given  by 

u  =  constant 

v  =  0 

r 


The  Rankine  vortex  is  obtained  by  replacing  the  central  part  of  a 
potential  vortex  with  a  core  which  rotates  as  a  solid  body.  Its  velocity 
profiles  are  given  by 

u  =  constant 
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2n62 


r  <  6 


w 


r 

2  nr 


r  >  6 


The  Rankine  vortex  has  been  used  by  Jones  [34]  and  Keller,  Egli,  and  Exley 
[41],  among  others,  for  vortex  stability  analysis  and  bursting  prediction. 

Lamb's  vortex  was  derived  as  a  solution  for  the  decay  of  what  is 
initially  an  infinite  line  vortex.  The  time  dependent  velocity  profiles  are 
given  by 


u  =  constant 


v  =  0 


w  = 


r 


2nr\ 


1  -  e 


4W 


Hall's  vortex  was  derived  as  a  solution  to  the  Euler  equations  for  a 
slender  conical  vortex  core.  The  velocity  profiles  are 


a  =  C(  K-\n- 
6 


90 


V 


--C- 
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w  -  A"  +  ^  -  In  ^  j 

where  K  is  a  shape  factor  for  the  profile  and  C  determines  the  magnitudes  of 
the  velocities.  This  model  was  used  by  Hall  [25],  Anders  [3]  and  others  for 
bursting  studies.  The  singularity  at  the  center  must  be  avoided  with  a 
patched  solution  to  a  diffusive  subcore. 

Mager  defined  his  velocity  profiles  as  polynomials  which  matched  values 
and  slopes  with  a  potential  external  flow  at  the  core  edge.  Only  the  u  and  w 
profiles  were  defined.  The  profiles  for  r  <6  are 


For  r  >  6,  the  Mager  vortex  velocities  are  the  same  as  for  a  potential  vortex. 

Several  modifications  of  the  Burgers  vortex  have  been  used.  The  Q 
vortex,  first  proposed  by  Lessen  [36],  uses  the  Burgers  vortex  circumferential 
velocity  profile  but  a  non-uniform  axial  profile.  These  profiles  are 

rj 

u  -  u0  e 
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where 


Q 


f 

2  n6u0' 


Several  other  investigators  have  used  this  vortex  model.  In  addition,  Powell 
and  Murman  [29]  derived  a  conical  analog  to  the  Burgers  vortex  which  is 
represented  as  a  system  of  ordinary  differential  equations. 


APPENDIX  B 


IMAGE  VORTEX  VELOCITY  DERIVATION 


Consider  first  the  two-dimensional  case.  The  primary  vortex  and  three 
image  vortices  are  assumed  to  be  infinite  and  parallel.  The  coordinate 
systems  and  definitions  used  in  Chapter  3  apply.  The  pertinant  geometry  for 
this  case  is  shown  in  Figure  32.  Note  that  for  simplicity  in  defining  angles, 
this  figure  views  the  crossflow  plane  from  the  upstream  direction,  while 
Figures  8,9,10  and  11,  for  instance,  view  the  problem  from  downstream.  For 
these  assumptions,  the  stream  function  in  the  XYZ  coordinate  system  due  to  the 
image  vortices  is  given  by 


(rsin0  +  2/r)2  =  (r2sin02)2 
(rcos0)2  =  (r2cos02)2 
(rsin0)2  =  (r3sin03)2 
(rcos0-2Zz)2  =  (r3cos03)2 
(r  sin  0  +  2 )2  =  (r4sin04)2 
(rcos0-2Zz)2  =  (r4cos04)2 

These  equations  can  be  solved  for  r2,  r3,  and  r4  in  terms  of  r  to  obtain 
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Figure  36.  Inage  vortex  geometry. 


r2  =  (r 2  +  4  lYr  sin  0  +  4  l2  )2 

l 

r3  =  [r2  -  4 1  zr  cost)  +  4  l2z)2 

I 

r4  =  (r 2  +  4 /Kr  sin  0  -  4  lz  cos  0  +  4 +  4 1]  )2 


The  stream  function  for  the  image  vortices  can  therefore  be  written 


¥,  =  —  (  - -In  (r  2  +  4  lYr  sin  0  +  4/1)  -  ^ln(r2  -  4lz  rcos  0  +  4  / 


2n\  2 


+  —  Inf  r2  +  4  lyr  sinQ  -  4 1  zr  cosd  +  4  ll  +  4  ll) 
4 


Then  the  induced  velocities  in  a  cylindrical  coordinate  system  fixed  at  the 
present  location  of  the  vortex  center  would  be 


\dYt 
V'  r  39 


r 


2  Zy  cos0 


2  lz sin  6 


2n\r2  +  4  lYrs\v\9+  4  ll  r2  -  4  lzr  cos0  +  4/1 


2 lY  cos  9  +  2  1 7 sin  9 


-L( _ 

2n\r2  +  4  lYr  sin  0  -  4  lzr  cos  9  +  4ly  +  4  lz 


w,  = 


dr 
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r  +  2/Ksin  9 


r  -  2/7cos0 


2n\r2  +  4lYrsin9  +  4ly  r 2  -  41  zr  cos  9  +  4l\ 
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2n 
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r  +  2  / y  sin  0  -  2  /z  cos  0 


\r2  +  4/  yrsin0-  4/zr  cos  9  +  4  lY  +  4 1  z 
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But  the  xr0  coordinate  system  is  moving  with  the  vortex  core,  so  the 
velocities  induced  at  r  -  0  must  be  subtracted  off  to  yield  the  image  vortex 
velocities  in  the  moving  coordinate  system  as 


f  2lycos9 

COS07  / 

2Zzsin0 

sin07 

\r2  +  4lYr  sin0  +  4l2 

2  ly  J  +  l 

^  r 2  -  4Zzr  cos0  +  4Z| 

2Zz  J 

r 


2  /y  cos0  +  2  Zz  sin  0 


ly  co s0  +  l,  sin  0\ 


2n\r2  +  4lrr  sin  0  -  4lzr  cos  0  +  4 12  +  4 1\ 


) 


u  im 


r  +  2/vsin0 


sinflA  f  r  -  2lzcos$  cos0^ 
2rt  \r2  +  4/Krsin0  +  4/2  2  lY  J  \r2  -  4 1  zr  cos9  +  4  lz  2lz  J 


r 


r  +  2/ysin0-  2  Zzcos0 


lY  sin  0  -  1 7  cos  0 


2n\r2  +  4  lYr  sin0  -  4Zzr  cos  0  +  4/2  4  l2z 


a 


Now,  in  the  three  dimensional  case,  the  image  vortices  are  not,  in 
general,  parallel  to  the  primary  vortex,  nor  to  each  other.  In  this  case,  the 
induced  velocities  must  be  multiplied  by  the  cosines  of  the  angles  between  the 
primary  vortex  and  the  image  vortices  to  get: 

r  (  2lrcos0  cos©^  f  2/zsin0  sin0v\ 

Vi~  2n^C°S<l> \r2  +  4  lYr  sind  +  4  l2y~  2  lY  j  +  C°S^r2  -  4/zrcos0  +  4l\ ~  2  lz 


2Zycos0  +  2Zzsin0 


r  ( 

+  — COS0yCOS0z  — z —  ,  , 

2  n  \r2  +  4Zyr  sin  0- 4Zzrcos0  +  4Z2  +  4Zz 


l  y  cos  0  lzs\n9\ 


J 


r  ,  f  r  +  2lyS\nd  sin 0 
W,“  2n^COS^Yyr2  +  4  lYrsin9  +  41$  2  lY  J 
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cos^Zl  _2 
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r  f  r  +  2  lY  sin  0  -  2  Zzcos0 

—  COS^yCOS0z  — - - - - 
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APPENDIX  C 


WING  GEOMETRY  DATA 


Three  different  types  of  wings  are  used  in  the  present  work.  The  first 
group  are  delta  wings  with  leading  edge  sweep  angles  ranging  from  55  to  76 
degrees.  The  second  group  of  wings  are  highly  swept  and  tapered.  These  wings 
are  often  referred  to  as  arrow  wings.  The  third  group  are  two  cropped  delta 
wings  with  identical  planforms  but  different  airfoil  sections. 

Four  different  delta  wings  are  used.  All  have  very  thin  flat  plate 
airfoils.  Pertinate  data  for  these  wings  are  listed  in  Table  1. 


Sweep  Angle  Aspect  Ratio  Half  Span/Root  Chord 


Wing  1  55  2.800 

Wing  2  65  1.865 

Wing  3  75  1.072 

Wing  4  76  1.000 


Table  1.  Delta  wing  geometry  data. 


.7002 

.4663 

.2679 

.2500 


The  arrow  wings  had  leading  edge  sweep  angles  of  71.2  degrees  and  taper 
ratios  of  0.10.  Half  span  divided  by  root  chord  was  0.4545.  The  aspect  ratio 
was  1.653.  Arrow  wing  number  1  had  an  NACa  0012  airfoil.  Coordinates  for  the 
airfoil  section  of  the  second  arrow  wing  are  listed  in  Table  2. 


X 

Y 

0.0000 

0.0000 

0.1250 

0.3359 

0.2500 

0.4506 

0.5000 

0.6064 

0.7500 

0.7247 

1.0000 

0.8182 

1.5000 

0.9520 

2.5000 

1.1191 

5.0000 

1.3448 

8.5000 

1.4809 

10.0000 

1.5195 

12.5000 

1 . 5444 

15.0000 

1.5630 

17.5000 

1.5720 

20.0000 

1.5813 

30.0000 

1.6214 

40.0000 

1.6398 

45.0000 

1.6282 

50.0000 

1.5901 

60.0000 

1 . 4344 

65.0000 

1.3121 

70.0000 

1.1627 

72.5000 

1.0792 

75.0000 

0.9921 

77.5000 

0.9006 

80.0000 

0.8069 

85.0000 

0.6132 

90.0000 

0.4156 

95.0000 

0.2153 

100.0000 

0.0113 

Table  2.  Arrow  wing  number  2  airfoil  coordinates. 

The  cropped  delta  wings  had  65  degree  swept  leading  edges  and  taper 
ratios  of  0.15.  Half  span  divided  by  root  chord  was  0.3964  and  aspect  ratio 
was  1.379.  One  wing  used  a  NACA  0012  airfoil.  The  other  used  a  NACA  65A005. 
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APPENDIX  D 


EVALUATION  OF  MODEL  CONSTANTS 


Consider  equation  (9c).  Three  different  w  velocity  profiles  are  chosen 
and  used  to  integrate  (9c)  to  obtain  a  relationship  between  pressure  at  the 
vortex  center  and  pressure  at  the  core  edge.  These  profiles  are  the  Rankine 
vortex,  Mager's  vortex,  and  the  Rankine  vortex  with  image  vortex  induced 
velocities  added  to  it. 

For  the  Rankine  vortex,  equation  (9c)  is  integrated  from  the  vortex 
center  to  the  vortex  edge  to  yield 


r2  r 6 

- r-—  r  dr 

4/r26  Jo 


or 


p(x ,6) - p(x ,0) 


— m 

4n262\2j 


rvn  _  r2 

8n\n62)  8na 


p(x , 0) 


P(x ,6)- 


r 2 

8  na 


so 


K 


a 


8  it 


and 


K, 


0 
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Next,  using  the  Mager  w  velocity  profile 


r 2  f6f  4  r  4r 3  <  rs\  7T2 

4/r262Jo  V<52  <54  66J  24/ra 


so 


p(x,0)  -  p(x  ,6) 


7  r2 

24/ra 


and 


=  0 


Finally,  for  the  Rankine  vortex  plus  image  vortex  velocities,  it  is  first 
assumed  that  within  the  vortex  core,  r  is  enough  smaller  than  1Y  and  lz  that  r2 
can  be  neglected  compared  with  41y2  and  41z2  in  equation  (10b).  It  is  further 
assumed  that  flattening  of  portions  of  the  vortex  core  which  come  close  to  a 
solid  wall  reduces  the  effect  on  core  pressure  of  the  velocities  induced  by 
the  image  vortex  on  the  other  side  of  that  wall.  This  occurs  because 
flattening  of  the  core  reduces  the  curvature  and  hence  the  radial  pressure 
gradient  of  that  portion  of  the  core  where  image  vortex  velocities  would 
otherwise  increase  the  radial  pressure  gradient.  As  a  result,  sine  and  cosine 
terms  in  (10b)  are  also  neglected,  and  the  expression  which  is  integrated  is 


where  C#  is  given  by  equation  (12).  The  integral  can  be  written 
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r 2 

4/r262 


(l-62C, 


r 2 

4/r262 


(  1  -26  2C,+64C2 


1 

2 


For  6  <  1Y  and  6  <  lz,  the  term  6*C#2  is  small  compared  with  62C# ,  so  the 
expression  for  pressure  is  approximated  by 


p(x ,0) 


p(x,6) 


r2 

8  na 


r 2 

4/ra 


62C 


* 


p(x,0) 


p(x,6) 


r2 

8/ra 


♦ 


a: 


8/r 


and 


1 

4/r2 


Next,  consider  equation  (9d) .  First,  for  a  Rankine  vortex,  v  «  0  and  the 
radial  variation  of  w  is  linear,  so,  if  u  is  not  zero,  (9d)  reduces  to 


dw 

Jx 


0 


which,  since 


rr  _  £r 
2nd2  "  a 
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can  only  be  satisfied  for  all  r  inside  the  core  and  all  x  if 


-  =  constant. 

aO) 

However,  equations  (9a)  and  (9b)  can  only  be  satisfied  for  this  case  if  vortex 
strength  and  core  area  are  both  constant  for  all  x.  A  more  useful  vortex 
model  can  be  obtained  by  using  the  Rankine  vortex  w  profiles,  but  allowing  v 
to  be  determined  by  integrating  equation  (9a)  with  respect  to  r.  The  u 
profile  will  still  be  modeled  as  uniform  in  the  core  but  not  necessarily  equal 
to  the  axial  velocity  outside  the  core.  Since  this  profile  is  axisymmetric , 
(9a)  can  be  rewritten 


d(ru)  t  d(rv) 
dx  dr 


so 


rv 


d(TlU) 

dx 


d  n 


where  n  is  a  dummy  variable  for  r.  Since  u  -  Uq  everywhere  in  the  core, 


rv 


du0( r2\ 
~dx\2  ) 


r  \du o 

2  J  dx 


Then,  in  equation  (9d)  the  Rankine  vortex  w  profile  still  causes  the  viscous 
term  to  vanish,  and  axisymmetry  makes  circumferential  derivatives  zero,  but 
the  other  terms  remain  to  give 


dw  f  r  \du.0duj  1  da0 

un - - - w 

0 dx  \2 )  dx  dr  2  dx 
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or 


/  r  dr 

rr( 

rr 

duQ 
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jx67\ 
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dx 

This  expression  is  integrated  over  the  area  of  the  core  to  obtain 


\_d£_  rI^>| 

2 dx  6 dx J 


rdu  o 
2~dx 


0 


1  da  1  dT  1  du0 

adx  r dx  u0  dx 


K,  -  1  ,  K2  -  1  and  K3  -  0 


If  the  Mager  profiles  are  used,  the  u  distribution  is  a  function  of  the 
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ratio  of  u  at  the  vortex  center  to  u  at  the  core  edge.  If  B  is  used  to 
represent  this  ratio,  which  Mager  calls  the  axial  velocity  form  parameter, 
then  u  is  given  by 


u 


“a 


B  +  (  1  -  B) 


and  the  v  obtained  by  integrating  (9a)  is 


duJrB  ,,  DJ3(rY  8i1 

V  dx  \  2  +r(  B\ 2UJ  6 . 


■  ru 


dx  \  2  2 


ra6(  1  -B)  —  {  3|  7 


dBj  l__3f  r 


When  equation  (9d)  is  integrated  over  the  area  of  the  core,  the  result  is 


\dT(2  B\  \  duJ\2\B-37\  f  1  i  \dB  I  da  f  157 1  B  -  227\  =  8  n 

fdx\ 5  15y  usdx\  630  )  \70jdx  2adx\  2520  )  3  Ree2au6 


which,  for  B-l  all  along  the  vortex,  reduces  to 


i  da  =  7  dT  1  du0  |  10 n 

adx  4T dx  u0  dx  Ree2au0 

K  \  ~  ,  ^2*1  and^3  =  lOn 
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